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ABSTRACT 



Four techniques for the numerical solution of partial 
differential equations and eigenvalue problems were 
Investigated. Typical problems considered were elliptic 
partial differential equations of the form 



U + U 
XX yy 



f (x,y) , 



( 1 ) 



or 

“yy 

where appropriate boundary conditions are specified so 
that the problem Is self-adjoint. 

The four methods are relaxation, Galerkln, Raylelgh- 
Rltz, and dynamic programming combined with Stodola's 
method, for eigenvalue problems. 

The results Indicated that for eigenvalue problems 

I 

relaxation or dynamic programming modified Is to be 
preferred usually and for partial differential equations 
Galerkln or dynamic programming Is preferred. 
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I . INTRODUCTION 



Initially for this thesis it was planned to Investigate 
methods for finding eigenfunctions and eigenvalues, with 
a particular interest in the oscillation in basins such as 
harbors and bays. The report by Angel [ 10 ] introduced 
a new method, that of dynamic programming, for the solu- 
tion of some partial differential equations. This method 
seemed promising and it was then extended here to the solu- 
tion of other partial differential equations. It also 
made it possible to Invert a linear differential operator 
and hence apply Stodola's method in finding eigenvalues 
and eigenfunctions. This led naturally to a comparison 
with several procedures already known to try to compare 
convergence, speed, and accuracy, by applying them to the 
solution of several rather simple problems. 

It is assumed that the reader has some familiarity with 
the techniques of replacing differential equations by dif- 
ference equations and the standard technique of separating 
variables in linear partial differential equations. The 
regions and their boundaries were assumed to be "nice," 
not excessively Irregular. Some knowledge of calculus 
of variations also is desirable but not necessary; ref- 
erence £5} provides more than adequate background. 

Pour problems were divided into two categories. In the 
first category were eigenvalue problems. Two typical ones 
were 
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U + U = - A^U 
XX yy 



subject to the constraint 



In A, 



( 1 . 1 ) 



U(x,y) = 0 



on 6A, 



( 1 . 2 ) 



where 6A is the boundary of the domain A; and second 



V“U = A^U 



subject to the constraint 



U(x,y) = — — ■■■ = 0 



in A, 



( 1 . 3 ) 



on 6A . 



(l.i+) 



In the second category were equations such as Poisson's 
equation 



U + U = f(x,y) 
XX yy ^ ' 



subject to the constraint 



in A, 



( 1 . 5 ) 



U(x,y) = g(x,y) 



and the biharmonic equation 



V“U = f(x,y) 



subject to the constraints 



on 6 A , 



in A, 



and 



U(x,y) = h(x,y) 



3 U f . 

= q(x,y) 

9n2 



on 6 A 



( 1 . 6 ) 



( 1 . 7 ) 



( 1 . 8 ) 



on 6A . 



( 1 . 9 ) 



Problems in the first category were numerically solved 
by a relaxation method, the Rayleigh-Ritz method, and 
dynamic programming combined with Stodola's method. 
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Problems in the second category were solved by Galerkin's 
method, the Rayleigh-Rltz method and dynamic programming. 

The domain used throughout this paper for comparisons 
was the unit square, with the constraint 

U(x,y) =0 on 6 A. (1.10) 

Computations were also carried out for L-shaped and trian- 
gular regions and for other boundary conditions. 

In the relaxation method and dynamic programming method 
in solving eigenvalue problems the initial estimates to 
the eigenfunctions had a special property. The first ap- 
proximation to Ui was picked so that it had no negative 
value in the domain. The approximation to U2 was picked 
so that it had negative and positive values in the domain 
and in addition it was orthogonal to Ui. The initial ap- 
proximations to U3 was picked in the same way as the one 
for U2 except that it had to be orthogonal to both Ui and 
U2. It may be difficult to make a suitable choice for some 
of the higher modes. 
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II. RELAXATION METHOD 



For many years relaxation techniques have been used 
to solve differential equations with and without the 
aid of computers. They are basically Iterative proce- 
dures In which a new approximation Is obtained from a 
previous approximation and Its residuals. 

A. DERIVATION OP EQUATIONS 

In this section a typical problem Is posed and solved 
by a relaxation method. 

Suppose the problem to be solved has the following 

form 

Z, , = Z + Z In A, (2.1) 

tt XX yy ’ 

where the unknown function Z must satisfy the differential 
equation In a simply connected region A In the xy plane, 
and for t>0. In addition the function Z Is required to 
vanish at points on the boundary 6A of the region A, 

Z(x,y,t) =0 on 6A. (2.2) 

A typical problem Is that of a vibrating membrane. The 
function Z(x,y,t) denotes the vertical displacement of 
the membrane. 

Now assume that the displacement has a representation 
of the form 

Z(x,y,t) = T(t) U(x,y). (2.3) 
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When Eq. (2.3) is combined with Eq. (2.1) the variables 
may be separated and a new equation is obtained of the 
form 



ijiii 

T“ 




+ U ) 

yy 



u 



-A- 



(2.4) 



This is equivalent to two equations 
T" + = 0 



(2.5) 



and 



U 

XX 



+ U 



yy 



-A^U, 



( 2 . 6 ) 



each with a parameter A. Further U must satisfy the 
boundary condition 

U(x,y) =0 on 6 A. (2.7) 

This is a typical eigenvalue and eigenfunction problem: 

to find the values A , or A , and the associated functions 

n 

U, or satisfying Eq. (2.6) and the constraint (2.7). 

The values A^, Ai ^ A 2 S A 3 < •••, are called the 
eigenvalues. With each eigenvalue is an eigenfunction 
U^. In this problem the eigenfunctions associated with 
different eigenvalues are orthogonal on the region A. 

Each eigenfunction is also called a mode. 

Consider a thin elastic membrane of a particular form, 
such as a very thin uniform sheet of rubber. Assume that 
the membrane is made fast at the boundary, while it is 
tightly stretched over the region with uniform tension. 

Also assume that damping is negligible. Then if an interior 
region of the membrane is pushed in a direction perpendicular 
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to the plane of equilibrium, it becomes distorted into a 
curved surface. The resulting area can be computed as 



S = ///I + dxdy, (2.8) 

A 

Assume that U and U are very small; Eq. (2.8) becomes 
X y 

approximately 

S = 11(1 + + igUy2) dxdy. (2.9) 

A 

The Increase in the area of the membrane due to the distor- 
tion is therefore approximated by 

6S = // (1 + 2 2 ) dxdy - // dxdy (2.10) 

X y 

A A 

= //(U^' + Uy 2 ) dxdy. 

A 

Hence the potential energy of the membrane in the deflec- 
ted position is 

PE = v/2 //(U^2 + (2.11) 

A 

where v is the tension, assumed to be constant over the 
region . 

Now consider any particular eigenfunction or mode. 

It follows from the solution of Eq. (2.4) that the deflec- 
tion is a periodic function of time and may be expressed 
in the form 

Z(x,y,t) = U(x,y) sin At, (2.12) 

except for a phase shift. Thus Eq. (2.11) can be rewritten 
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as 



PE = ^ // (U 2 + u 2 ) dxdy sln^ Xt . (2.13) 

X y 

A 

The maximum value of the potential energy Is 

^^max " //^V (2.14) 

A 

The kinetic energy of an element dm = pdxdy of the membrane 
is 

^pdxdy(U^)2 = Jspdxdy(U2x2 cos^xt), (2.15) 

where p denotes the mass per unit area of the membrane. 
Therefore, the kinetic energy of the vibrating system is 

KE = hX^p //U^dxdy cos^xt, (2.l6) 

A 

and the maximum value of the kinetic energy is 

KEmax = //U^dxdy. (2.17) 

A 

If it is assumed that the energy is constant then, the 
maximum values of the potential and kinetic energy are 
equal for individual modes, and therefore 

^SX^P //UMxdy = v/2 J/(U^" + U^") dxdy (2.l8) 

A A 



or. 



V J/(U^2 + u^2) (2.19) 

A 

p J/U^dxdy 
A 
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The first eigenfunction Ui is the function v/hlch minimizes 
this quotient and the first eigenvalue Is the corres- 

ponding value for the quotient. The next eigenfunction 
is the function U 2 which minimizes the quotient in the 
space of functions orthogonal to Ui, and is the corres- 

ponding value of the quotient. The third eigenfunction 
U 3 minimizes the quotient in the space of functions ortho- 
gonal to Ui and U 2 , etc. The set Ui, U 2 , ... is unique 
except that if two or more eigenfunctions have the same 
eigenvalue, they may be replaced by linear combinations 
of themselves. It is seen that Xi < X 2 < X 3 ... 

If the mode U is known, Eq. (2.19) can be used to obtain 
X^. An Interesting fact is that a rather poor approximation 
to the first mode Ui, chosen to satisfy the appropriate 
boundary condition, will yield a surprisingly good estimate 
for Xi^. This result is apparently due to Lord Rayleigh, 
and the quotient is sometimes called Rayleigh's quotient. 



B. COMPUTATIONAL ROUTINE 

In this section the results of section A will be used 
to develope a method to solve Eq. (2.6). 

If Eq. (2.6) is expressed as a difference equation 
then it may be written as 



'“l.l.j " "1-1, J - 



h^ 



k‘ 



= “i.j 



( 2 . 20 ) 
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where h and k denote the x and y mesh size respectively. If 
the mesh sizes are equal then Eq. (2.20) can be rewritten as 



U. 



^ “l.J.l ^ “l-l.J ^ '2.21) 






^ - X^h^ 



The procedure used to solve the partial differential 
equation was to pick a function U° which satisfied the given 
boundary conditions as a first approximation to U(x,y). 

This function need not be a very good approximation to U, 
and in fact step functions were sometimes used. 

From Eq. (2.19) an approximation to X^ was obtained 
by assuming that the approximation UO was the desired func- 
tion U. With this first approximation to X^ Eq.(2.21) 
was used to obtain an improved estimate of U. The form 
of Eq. (2.21) used to obtain the vth approximation was 

v-1 V-1 V V 



(U.^, . + U. + U. T . + U. . t) (2.22) 

V 1+1, J i,J+l 1-1,0 1,0-1 

U = » 

(i| - X^h^ ) 



in which the y-lst estimate of X was used. By alternating 
Eq. (2.19) and Eq. (2.22) a close approximation to X^ and 
U were obtained. The solution converged to the smallest 
eigenvalue Xi and the associated eigenvector Ui . 



C. HIGHER ORDER EIGENVALUES 

In this section the method is extended to find the 
larger eigenvalues and eigenfunctions. 

To obtain the larger eigenvalues and eigenvectors the 
same procedure as that in section B was used except that 
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additional equations were added to force the eigenfunctions 
to be orthogonal to those already found. Define for 
1 = 1 , 2 , . . . , n to be the 1^ eigenfunction, associated 
with All higher eigenfunctions must be orthogonal to 

the lower ones obtained. Orthogonalizatlon was accomplished 
by subtracting out multiples of the lower eigenfunctions 
already obtained. This was effected by expressing U as 



^ 1+1 ~ ^ 1+1 ^ * 
^ new ^ -^old j=l ^ ^ 



( 2 . 23 ) 



where 



//Ui+i U dxdy 



A old 



C. = 



, J=l,2,...,i. (2.2i|) 



^ dxdy 



For each eigenfunction desired a different initial 
approximating function was used. The method gave the eigen- 
values and associated eigenfunctions in numerically increas- 
ing order. The method was very simple and effective for 
lower eigenvalues. 

The fault of the method is that convergence was poor 
and computation times were large if the number of Intervals 
in both directions was large. 
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III. DYNAMIC PROGRAMMING [8,10j 

The method of dynamic programming has recently been 
applied to the solution of partial differential equations, 
[10]. The difference equation may be regarded as leading 
to a system of linear equations with a large number of 
unknowns. The method effectively reduces the number of 
unknowns involved so that a number of systems are to be 
solved, each one of much lovjer dimension. In one case, 

for example, using a grid with N+1 intervals in the region 
in each direction, N systems each with N unknowns are solved, 
rather than one system with N - squared unknowns. 

In section A, the method is applied to the solution of 
Poisson's equation over a rectangle, following the paper of 
Angel ClO] • section B, the method is extended to a re- 

lated eigenvalue problem by combining it with Stodola's 
method. In section C, the method is extended to the solu- 
tion of the blharmonic equation, by applying dynamic pro- 
gramming twice. Finally, in section D, the method is 
applied to the eigenvalue program involving V^U = X^U, again 
by combining the method with Stodola's method. 

A. DYNAMIC PROGRAMMING FOR PARTIAL DIFFERENTIAL EQUATIONS 

Consider the solution of Poisson's equation 

U + U = f(x,y) in A, (3.1) 

XX yy > 

where U = U(x,y) is subject to given boundary conditions 
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such as 



U(x,y) = g(x,y) 



on 6 A. 



(3.2) 



It may be noted that Eq. (3.1) is the Euler equation 
associated with the variational problem 

I(U) = Min //(U^^ + U/ + 2fU) dxdy (3.3) 

U A ^ y 



where the function U is chosen from the class of functions 
with first partial derivatives belonging to over A, and 
satisfying the boundary conditions of (3.2) on <5 A. 

Let the region A be discretized by choosing n+1 and 
m+1 equally spaced points in the independent variables 
X and y respectively. Then Eq. (3.3) niay be rewritten in 
a discrete version with equal intervals as 



n m ( 



I(U) = Min I I 



“i,j 



(U. . - U. . T )2 



(3.4) 



+ (U. . - U. , .)2 + 2f. . U, .h^ 

i.J 1 - 1 , J i,J ij J 



where (u (u. ) {u and {U. } are determined 

^ SO i,m 

from the boundary conditions of Eq. (3.2). Now, if all terms 
in Eq . (3.4) involving only boundary values were removed, 
while not affecting the solution, a more convenient form 
of Eq. (3.4) is obtained 



I(U) = 



n 

Min I 

U. . 1=1 

i,J 



m 

I oil 1 - 

J=1 






(3.5) 



m-1 m-1 

+ I 2h“f ,“l , + I (U, , - U, , ,) 

1=1 ^ SO 1=1 ^ s^ sJ 
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In vector notation Eq. (3*5) may be rewritten in the form 



n _ _ _ 

I(U) = Min I (<QUp,Uj^> + + Sj^ (3.6) 

Ur 1=1 

+ <Ur - Ur.i, Ur - Vl > < ^^R>Ur» 

— T 

In this, Ur = (Ur , . . . , Ur ) . This relation now defines 

m-1 

a symmetric matrix Q, vectors Fr, and fR, and a scalar 

Sr by 



Q = {q. where q. . 



(2 1 = j 

-1 |i-j| = 1 
0 otherwise 



(3.7) 



•■r = '■r.I 



-2“r.o J = 1 
-2“R,m J = 

0 otherwise 



fR = {2h^fR^j} 

"" ^R,o ^R,m 



The notation as in <rR, Ur> , denotes the scalar product of 
the two corresponding vectors. The matrix Q remains constant 
while rR and Sr are functions of the boundary conditions 
only. 

Now in order to solve Eq. (3.6) a sequence of dynamic 
programming problems was considered. Let 



_ n _ _ _ 

F„ (V) = Min_ I (<QU, , U.> + <r, , U . > (3.8) 

«R---Vh=« 



+ s. 



. + <U. - U. U. - U. T> + <f., U.>) 

1 1 1-1’ 1 1-1 l’ 1 
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where U^, i = V' and O' is given by Eq. (3.2). Now Eq. (3-8) 
n— 1 n 

can be rewritten as 

Fj^(V) = Min V ^3.9) 



+ <U^ _ v,U^ - V> + <f^ , U^> + <QU^^., U^^3_> 



+ ... + s + <v - u , , tr - u ,> + <f , u >) 

n n n-1 ’ n n-1 n’ n 



This may be done since the minimization over 

U - can be commuted with the minimization over Urj. This 
n-1 R 

is a common technique of dynamic programming based upon the 
principle of optimality, [e], This allows Eq. (3-9) to be 
rewritten as 



Pj,(V) = Min(<QUp. Ujj> + <rp, U^> + Sj, (3.10) 

^R 



+ <Upj V, Up V> + <fp, Up> + Fp4.T (Up) 



R 



R’ R R+1" R- 



The final equation is 



F (V) = <QU ,U>+<r,U>+S 
n ^ n’ n n* n n 



(3.11) 



+ <U - V, U - V> + <f , u >, 
n * n n’ n ’ 



and is known from the boundary conditions Eq. (3.2). 
Since Fj^(V) is quadratic in V Eq. (3.11) may be rewritten 
in the form 



Fr(V) = <Aj^7, V> + <b^Y> + (3.12) 

Now by substituting from Eq. (3.12) into Eq. (3.10) and 
then differentiating with respect to an expression for 



19 



is obtained of the form 



= (I + Q + 



-1 



V - 



^R+1 ■*■ 



. (3.13) 



This is obtained by substituting related Eq. (3.13) into 
Eq. (3.10) and then combining Eq. (3.10) with Eq. (3.12), 
the various quantities in Eq. (3.13) are defined by these 
steps as 

Apj = T - (T + Q 



R 



'R+1 



^R ^ ^R+1 ^R " 



^R+1^ 


(3.14) 


7 ^^R-1 ^R ^R^ 


(3.15) 


(I t Q + Aj^^^)-1 


( 3 . 16 ) 


R+1 ^R . 





2 2 

with initial values determined from Eq. (3.11) as 



A = I, 
n ’ 



= -2“n- 



C -<(I + Q)U,U> + <r - U> + S. 
n n ’ n n n n 



(3.17) 

( 3 . 18 ) 

(3.19) 



The matrix (T + Q + A^^q^) is nonsingular, 1 . Thus it has 
inverses which may be computed beforehand, since it is de- 
pendent only upon the type of operator. 

Due to the fact that only the values of tl_ are desired 

n 

the quantities need not be calculated. 

The procedure is to calculate the quantities in 
Eq. (3.14) repeatedly until A 2 and ba are obtained. Then 
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U is found from Eq. (3 .13) with V = U^, which is known 
from the boundary conditions. Next, Eq. (3.13) is repeated- 
ly solved using the stored values of and bj^, and the 
last value of as V. 



Thus, the problem was solved by n-1 Inversions of sym- 
metric matrices of order m-2. While these matrices may be 
large there are efficient computer routines available to 
determine the Inverses. Once the Inverses are found, they 
may be stored for future use, since they are based only on 
the geometry of the region A. Thus for several problems 
over the same region the Inverses may be entered into the 
program as data. Also they have the property that if less 
than n+1 grid points are required a reduced number of the 
matrices may be used. 

B. DYNAMIC PROGRAMMING FOR EIGENVALUE PROBLEMS 

In the first section of this chapter it was seen that 
dynamic programming could be used to solve partial differ- 
ential equations. In this section, the program is modified 
to solve a related eigenvalue problem by Stodola's method. 

Consider the problem of the vibrating membrane consi- 
dered in Chapter II. The differential equation for the 
eigenvalues and eigenfunctions is 



where U = U(x,y) is subject to the boundary condition 



R 




in A» 



( 3 . 20 ) 



U(x,y) = 0 



on 6A. 



( 3 . 21 ) 
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This problem was solved by two related methods. In 
both Eq. (3.20) is regarded as a special case of Eq. (3.1) 



in which 



f(x,y) = -^^U(x,y) 



in A 



( 3 . 22 ) 



and 



f(x,y) = 0 



on 6A 



(3.23) 



1 . First Routine 

The difficulty here was that the function f was 
known only on the boundary and \ was unknown. The first 
step, to overcome this obstacle, was to choose a function 
U°{x,y)which satisfied the given boundary conditions. Then 
an estimate of was obtained, say (x°)^, by using Ray- 
leigh's formula. Next a new estimate of U, say U^(x,y), was 
obtained using dynamic programming to solve the equation 



By repeating this sequence of Rayleigh's formula and dyna- 
mic programming, a good approximation to the minimum eigen- 
value and the associated eigenfunction was obtained. The 
next two eigenvalues and vectors were also obtained by the 
process of forcing the higher eigenfunctions to be orthogo- 
nal to ones already obtained as was done in relaxation. 

The inverse matrices used in the routine were calculated 
in determining the first eigenvalue and eigenfunction; 
they were then stored so that it was not necessary to recal- 
culate them for subsequent eigenvalues. 




(3.24) 
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2. Stodola's Method 



The second form is more like the usual form of Sto- 
dola's method and hence a brief of Stodola's technique is 
given first. 

An initial function Vq satisfying the boundary con- 
ditions is selected. It may be considered to be of the 
following form 



Vo = aiUi + aaU 2 +. . . (3.25) 

where ai is not equal to zero. For convenience was 
normalized; the Loo norm of Vq, | |Vq| | , is defined as 

MAX |Vo(x,y)| = I |Vj | ; (3.26) 

A 

and I [VqI I is set equal to one. Now consider 



L 






(3.27) 



and 





m 




r ^ 1 


2m 






aiUi + az 


^ 1 
X 2 


U 2 + . . . 


V > 
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In Eq. (3.28) it can be seen that the relative size of the 
components Ug , U3,... are decreasing by a factor (Xi/Xa)^} 
(Xi/Xs)^*... respectively. After a few applications of the 
operator, L~\ to V^, the leading term will dominate. 

Thus the functions obtained will approximate Vj, except 
for a constant factor, and the ratio of successive iterates 
will approach a constant, 
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3 . Second Routine 



This was applied as follows . Let Z be defined by 

m 

the equation LZ„ = V n , so that 
^ m m-1 ’ 



Z 



m 



L 




(3.29) 



This equation was solved by dynamic programming. It was 
found convenient to normalize each Iteration. Let 



V = 
m 




(3.30) 



Then the approximation can be made 




(3.31) 



The sequence of functions Vq, Vi,... converges to Ui. 

The process was term.lnated when successive approx- 
imations were sufficiently close together, say, 

I Iv^(x.y) - |<6. (3.32) 



where e Is a preassigned small positive number. The result- 
ing function Is an approximation to Ui, and Xi Is ob- 
tained from Eq. (3.31). 



C. DYNAMIC PROGRAMMING FOR A HIGHER-ORDER OPERATOR 

In section A and B of this chapter It was shown how 
dynamic programming could be used to solve Poisson's equa- 
tion and the vibrating membrane program, respectively. 

In this section, by another modification to the routine, 
the blharmonlc equation can be solved. 



Assume the problem to be solved Is of the form 



V^CU) = f(x,y) 



in A, 



(3.33) 



where U = U(x,y) is subject to the constraints 

U = g(x,y) on 6A (3.3^) 

■ 

3^U/9n^ = P(x,y) on 6A. 

Clearly Eq. (3 *33) may be rewritten in the form 

$(x,y) = f(x,y), (3.35) 

where 

U(x,y) = $(x,y). (3.36) 

Now, since U and 3^U/3n^ are both known as the boundary, 
$(x,y) may be approximated on the boundary. On a rectangle, 
for example, the following relations determine $ on the 
boundary 

(3.37) 



$(0,y) = -P(0,y) 


+ 


^yy 


(0,y) 


•J>(n,y) = P(n,y) 


+ 


n 

yy 


(n,y) 


$(x,0) = -P(x,0) 


+ 


u 

XX 


(x ,0) 


^ $(x,m) = P(x,m) 


+ 


U 

XX 


(x ,m) . 



The usual finite difference scheme may be used to approxi- 
mate U and U and thus the relations of Eq. (3.37) may 
XX yy 

be approximated by 



(Un -j.-, - 2U„ . + U„ . T ) 

^ = _p^ ^ 2 - >J - J -— (3.38) 



0,j 
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(3.38) 




<I> 



i,0 



-P 



1,0 





P 



1 ,m 




for 1=1 



n-1, 0=1 



m-1. 



and '**Q Qj 0 ’ ^0 m’ m known from the boundary con- 



ditions. Thus <J>(x,y) Is now approximated on the boundary. 
First Eq. (3.35) Is solved by dynamic programming to obtain 
the function $ In the region. Then Eq. (3.36) Is solved 
by dynamic programming to obtain the desired solution. 

D. DYNAMIC PROGRAMMING FOR HIGHER ORDER EIGENVALUE PROBLEMS 
This method may be extended to the corresponding eigen- 
value and eigenfunction problem much as before. One such 
physical problem Is that of a vibrating uniform plate with 
hinged edges. 

Assume that the differential equation has the form 



The technique of sections B and C may be applied In two 
steps to the solution of the problem. Let $ be defined as 



V“U = 



In A 



(3.39) 



and the boundary conditions are 




on 6A 



(3.^0) 



$ = v^u 



(3.^1) 
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Then the two equations to be solved are 

(3.42) 

and 

y2yn+l ^ ^n^ (3.43) 

subject to the boundary conditions, Eq. (3.40). 

In order to start the routine an initial estimate for 

the function U(x,y), say {U..}, is made. The value of A 

^ J 

is estimated as was done in section B. Then Eq. (3.42) 

and Eq. (3.43) are solved by dynamic programming to get the 

next approximation {U. The same criterion for stopping 

1 > 0 

is used as that in section B. 
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IV. RAYLEIGH-RITZ METHOD [^,5,6,9] 



The Raylelgh-Rltz method has been used for many years 
to obtain approximations to the solution of partial differ- 
ential equations and eigenfunction problems. In this method 
the problem is posed as a minimization problem, say invol- 
ving an integral. Then some linearly independent functions 
which satisfy the boundary conditions are chosen. The solu- 
tion is approximated by a linear combination of these. 
Finally the coefficients in the approximation are chosen 
so as to effect minimization. This leads to an eigenvalue 
problem Involving symmetric matrices. 

The functions chosen may be, for example, polynomials 
of low degree, or trigonometric functions. Assume that 
homogeneous boundary conditions are given. Let = $j^(x,y) 
be n functions which satisfy these and approximate the solu- 
tion U by 

U(x,y) = I C . (^.1) 

k=l ^ ^ 

A. PARTIAL DIFFERENTIAL EQUATIONS 
In this section the solution of 

U +U =f(x,y) in A, (1.3) 

XX yy * 

is again considered. It is the Euler equation associated 
with minimizing the Integral 
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(4.2) 



//(U ^ + U ^ + 2fU) dxdy . 

A X y 

Substituting from Eq. (4.1) into Eq. (4.2) and integration 
results in a function which may be written 

I = I(C^,..., C^), (4.3) 

for functions of the form (4.1). The minimization of this 
function and an approximate solution to Eq. (1.3) is thus 
obtained by solving 

= 0, k = 1, 2, . . . , n. (4.4) 

The effectiveness of the procedure of course depends upon 
the choice of the approximating functions $j^(x,y). 

B. EIGENVALUE PROBLEMS 

The method is also applicable to eigenvalue problems. 
Consider again the problem in Chapter II of finding eigen- 
values and eigenfunctions for the equation 

U + U = -X^u in A (2.6) 

XX yy 

subject to the boundary condition 

U(x,y) =0 on 6A. (2.7) 

2 

In many problems Xi, the lowest eigenvalue, is the minimum 
of the ratio of two Integrals. This fact was shown in 
Chapter II and Eq. (2.19). 

If an approximation to U is chosen as in section A, since 
each function satisfies the linear homogeneous boundary 
condition, the sum shown in Eq. (4.1) satisfies it. If 
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this sum is substituted Into Eq. (2.19) the resulting values 
define an upper bound for X, for all choices of the C’s. 
Further also, the C's may be chosen to give a least upper 
bound over the subspace spanned by <J>i, $ 2 ,..., 

If that approximation for U is substituted into Eq.(2.19) 
the numerator and denominator become quadratic forms in 
(C C ) = U. The numerator has the form 

i ? a. C C = C^A C (4.5) 

i=l j=l ^ 

where 

a 4 . = + '^>. '^>. ) dxdy, (4.6) 

A ^x ^x ^y ‘Jy 

and the denominator has the form 

n n . 

I I b..C.C. = C (4.7) 

1=1 j=l J 

where 

b. . = //$. $ dxdy. (4.8) 

-^J A ^ 

These define the matrices A and B. The first eigenvalue 
Xi minimizes the quotient in Eq. (2.19) and hence the mini- 
mum value of 

C^A C/C^B C 

defines the minimum value of the quotient in the subspace 
s.panned by $ 1 ,..., 
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To find this is equivalent to minimizing the quadratic 



form 




X 1 ^ = Min A C 

C 


(4.10) 


subject to the constraint 
assumes the value one 


that the second quadratic form 


C^B C = 1. 


(4.11) 


The problem may be solved 


in two steps. First find the eigen 



values and eigenvectors of B. Let the eigenvalues of B 
be i=lj 2 ,..., n, and the associated normalized eigen- 

vectors V^. Let T be the matrix 

T= (Vi, V 2 ,..., V^) diag(l/wi,..., 1/w^). (4.12) 



Then the transformation 




C = T D 


(4.13) 


reduces the constraint (4 


.11) to the form 


^^B C = D^D = 1. 
condition (4.10) becomes 


(4.14) 


Xi^ = Min D^E D, 
D 


(4.15) 


where 

E = T^A T, 


(4.16) 



subject to the constraint (4.14). 

Hence the value of is the smallest eigenvalue of 

the matrix E. 
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Let Di be the associated normalized eigenvector of E. 
Then the corresponding minimizing function has coefficients 
determined from Eq. (4.13) 

Cl = T Di . (4.17) 

This value of Ai^ls an upper bound for the first eigen- 
value, the least upper bound In the space of functions 
spanned by i , . . . , . In a similar way the second eigen- 

value of E furnishes an upper bound to the second elgenfunc 
tlon, etc. 

For simple problems, at least. It seems easy to choose 
the 4>’s so that a good bound for the lowest eigenvalue re- 
sults. It Is not clear how to make good choices to get 
good estimates for the higher eigenvalues . 
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V. GALERKIN’S METHOD [^,5,6j 



Sometimes the solution to boundary-value problems In- 
volving partial differential equations can be obtained by 
forming an associated functional In the form of a definite 
Integral which Is to be made stationary. For a problem of 
this type, the Raylelgh-Rltz method Is often an effective 
procedure for the determination of an approximate solution. 
However, In many Instances It Is difficult to find this 
functional. In such cases, the Galerkin method is often 
effective . 

A. PARTIAL DIFFERENTIAL EQUATIONS 

Consider the linear homogeneous boundary value problem 
of the form 

LU(x,y) = f(x,y), (5.1) 

subject to linear homogeneous boundary conditions. The 
symbol L stands for a linear differential operator, such 
as • 

Suppose that an approximate solution is taken In the 

form 



U^(x,y) = C^$j^(x,y), (5.2) 

where the coefficients C^,..., are constants, to be de- 
termined, and, as in the previous chapter, the functions 
$ 1 ,..., are picked to satisfy the homogeneous boundary 
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conditions. The coefficients are dependent upon the number 
of functions picked and therefore, must be recomputed if 
a larger number of functions are later chosen. V/hlle the 
function satisfies the boundary conditions, it will not, 
in general, satisfy Eq. (5.1). Thus there is a residual, 

LUn(x,y) - f(x,y) = R^(x,y), (5.3) 

which can be viewed as an error or penalty function. Assume 
the solution for U(x,y) can be expressed by an infinite 
complete series of these linearly independent functions 
in the form 



U(x,y) = I c.<l'. (x,y). (5.^) 

k=l ^ ^ 

Now U^(x,y) in Eq. (5.2) represents a sequence of partial 
sums which approximate U(x,y). Now if the condition is 
imposed that L(U^) - f be orthogonal to each function 
'^’^(x,y) on the demain A, the following set of equations 
are obtained 

//(l(U^) - f) dxdy = 0, for k = 1, 2,..., n 

An K 

or by use of Eq. (5.3) (5.6) 

//r (x,y) <^^(x,y) dxdy = 0, for k = 1, 2 ,..., n. 

A ^ K 

If Eq. (5.6) is to hold as n “ it follows that 

11m = 0 (5.7) 

n-><» 



3 ^ 



Let R(x,y) be an arbitrary function satisfying the 
homogeneous boundary conditions. Since for k = 1, 2,... 
forms a complete set of functions^ constants a^,..., a^,... 
can be found such that 



R(x,y) = I a,$^^(x,y). 
k=l ^ ^ 



(5.8) 



Thus , 

// lim R„(x,y) n(x,y) dxdy = 0, (5.9) 

A n-»-“ ^ 



for any arbitrary R(x,y), so that by the fundamental lemma 
of variational calculus Eq. (5.7) is true "almost every- 
where". Suppose further that L(U^^) L(U)j then Eq. (5.1) 
holds, by the use of Eq. (5.7). 

Galerkin's method requires that the error function 
Rn(x,y) be orthogonal to each of the functions or that 

Eq. ( 5 . 5 ) and Eq. (5.6) hold. Now by substituting Eq . (5*2) 
into Eq. (5.5) an integral is obtained of the form 



// 

A 



n 



L( I 

k=l ^ ^ 



- f 



dxdy = 0, 



( 5 . 10 ) 



i = 1, 2,..., n. 

This results in a system of n linear algebraic equations 
in n unknowns C^,..., C^. Furthermore, the system is in- 
homogeneous unless the function f(x,y) is orthogonal to 
each $j^(x,y). Eq. (5.10) may be rewritten 
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(5.11) 



n 

I C //l($ ) $ dxdy = j/f4> dxdy, 

j=l ^ A J ^ A ^ 

k = 1, 2,..., n. 

In matrix notation this becomes 

B C = P, (5.12) 

where 

C = (C^ C 2 . . .C^)^ (5.13) 

B = {b^j} 1 = 1, 2, . . . , n, (5.14) 

<J “ 1> 2,..., n 

F = (f^,..., f^)^ (5.15) 

and 

b.. = //l($.) $. dxdy (5.16) 

-LJ A ^ 

f. = //f$. dxdy (5.17) 

^ A ^ 

The constants are obtained by solving the system Hq. 

Eq. (5.12). 

Collocation, a convenient variation. One way to get 

an approximate solution for the C's is the following. 

Choose n points of the region A. Evaluate the terms in 

the integral of Eq. (5.10) at these points. This yields n 

equations for the unknowns C . This method, called colloca- 

k 

tlon, is not generally so accurate but it is quicker than 
carrying out the integrations to define the coefficients 
Eq. ( 5 . 16 ) and Eq. (5.17) . 
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Galerkins method Is also useful to obtain a direct 



solution of variational problems. Assume it is desired to 
minimize a functional of the form 

I(V) = //P(x,y V,V^,V ) dxdy. (5.18) 

A ^ ^ 

It is desired to find an extreme value for the functional 
subject to the condition that V(x,y) is prescribed on the 
boundary 6A of the domain A. It is known that if the exist- 
ence of an extremizing function U(x,y) is assumed and that 
the function F possesses continuous derivatives of the 
second order with respect to its arguements, there results 
the condition 



//h(x,y) 

A 





T~ 



dxdy = 0 



(5.19) 



for an arbitrary function ri(x,y) which has piecewise con- 
tinuous derivatives in A and which vanishes on 6A. This 
is derived by considering a function of the form 



V(x,y) = U(x,y) + en(x,y). 



( 5 . 20 ) 



and differentiating with respect to e . If n(x,y) is other- 
wise arbitrary, then one form of the fundamental lemma of 
variational calculus requires that 



P 



U 





0 . 



( 5 . 21 ) 



Suppose now that n(x,y) is the kth function $j^(x,y) and 
that an orthogonality requirement is imposed as 
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Pu^ - PUy) dxdy = 0, (5.22) 
for k = 1, 2,..., n. 

Plnally, assume that the function U(x,y) which effects the 
minimization can be represented satisfactorily by a finite 
series 



Uj^(x,y) = I Cj^$j^(x,y). (5.23) 

k“l 

Eq. (5.22) then defines a system of n algebraic equation 
to be solved for the n unknowns c^,... c^. Thus Galerkin's 
method is applicable in the solution of variational problems. 
However, much of its value lies in the fact that it is not 
necessarily connected with a variational procedure. 



B . EIGENVALUE PROBLEMS 

Suppose that it is desired to solve an eigenvalue prob- 
lem of the form 

LU = XU(x,y) in A, (5.2^^) 

and 



U = 0 



on 6A. (5.25) 



Assume as in the first section that the function U(x,y) 
can be approximated in the form of Eq. (5.23). The problem 
is solved by considering equations of the form 



// 



A 



L(U ) - XU 
^ n"^ n 



$ . dxdy = 0 



( 5 . 26 ) 



for j = 1, 2,..., n. 
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Now by substituting for the functions a system is ob- 
tained of the form 



n 

I 



j=l 






0 



( 5 . 27 ) 



for i = 1 , 2 ,..., n. 



where 



a 



i J 




( 5 . 28 ) 



n,j = 



( 5 . 29 ) 



This has nontrivial solutions for the C’s if, and only if, 

' 5 . 30 ) 

where A is the determinant of the matrix. Thus the values 
of may be found by solving the characteristic equation 
(5.30). For each eigenvalue there corresponds a system 
of equations ( 5 . 28 ) to be solved for the eigenvector 



n 

I 

k=l 



I (a 



k,j 









( 5 . 31 ) 



for j = 1 , 2 ,..., n. 

Now for each eigenvalue X^ this system of n homogeneous 
equations may be solved to give the values where 

this coefficient corresponding to the i^ eigenvalue. 

Thus the eigenvalues are obtained, with their correspond- 
ing eigenvectors. The functions should be chosen to 
reflect whatever characteristics the solution is felt to 
have . 
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VI. DISCUSSION OF VARIOUS METHODS 



FOR SOLVING EIGENVALUE PROBLEMS 

In Chapters II, III and IV, three methods were developed 
for finding the first three eigenvalues and eigenfunctions. 
All three gave satisfactory values for the eigenvalues and 
eigenfunctions, and satisfactory times for the rather simple 
problems considered here. These methods were used to obtain 
solutions of the following two problems: 

U + U = - X^U in A, (6.1) 

subject to the constraint 

U(x,y) =0 on 6A; (6.2) 

and 

V‘*U = X^U in A, (6.3) 

subject to the constraints 

U(x,y) = ^ =0 on 5A. (6.4) 

9n^ 

In this chapter the computation and numerical results are 
discussed and compared. 

In all of the methods a set of functions was needed. 

In the Raylelgh-Rit z method these were the basis for the 
approximating functions; in the other two methods they 
were the initial estimates of the functions. The ones usu- 
ally chosen were 
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(6.5) 



Ui° = (x - x^) (y - y^) , 

Ua° = (x - x^)(y - - x), (6.6) 

U3° = (X - x")(y - - y). ’ (6.7) 

Step functions were also used as first estimates in the 

iterative methods; these increased the number of iterations 
required some, but not much, particularly for the higher 
modes . 

Usually the range of the Independent variable was di- 
vided into ten equal sub-intervals, in going to a differ- 
ence equation. In the Rayleigh-Rit z method this did not 
yield sufficient accuracy and it was found necessary to go 
to forty intervals. Twenty-five intervals were also used 
in some computations; the intermediate number was chosen 
because of storage requirements in dynamic programming. 

In the iterative procedures some convergence or stop 
criterion was needed. When a relaxation procedure was used 
together with Rayleigh's formula for estimating the eigen- 
value, computation was terminated whenever the eigenvalue 
did not decrease by at least 0.002. In Stodola's method, 
the routine was terminated whenever the norm of the change 
in the function U(x,y) was less than 0.002. 

The relaxation method required the least time for this 
simple problem. Most of the time required for dynamic pro- 
gramming was spent in inverting the matrices. Dynamic pro- 
gramming yielded a very accurate approximation to the 
eigenfunction . 
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For some reason it was necessary to use forty Intervals 
to get the desired accuracy with the Rayleigh-Ritz method. 
When ten intervals were used, the eigenvalues of the matrix 
E were significantly too small, apparently due to errors 
in the integration and transformation routines. There were 
at least two other disadvantages of the Rayleigh-Ritz method. 
First it was tedious to program. Second, there may be some 
difficulties in choosing the functions particularly 

if some of the higher modes are desired. The results of a 
set of computations are given in Computer Output 5 , in- 
cluding the three eigenvalues, the associated coefficients 
and the values of the corresponding functions at various 
points. The routine is shown in Computer Program 5. The 
necessary matrix transformations and solutions for the 
eigenvalues were carried out using programs TRED2 and TGL2 
respectively, [?]• 

The simple relaxation method of Chapter II had the ad- 
vantage of being the simplest to program and to run. It 
had the disadvantages that terminal convergence was slower 
than in the method of dynamic programming and if a large 
number of points were Involved the computing time increased 
greatly. The results of a set of computations are given 
in Computer Output 1. The routine is shovm in Computer 
Program 1. 

Dynamic programming had the advantage of yielding very 
accurate values in a small number of iterations. The dis- 
advantages were that it was relatively difficult to program 
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and that It took quite a bit of time to generate the inver- 
ses of the matrices required. Computer Output 7 shows the 
results by this method and the program used Is in Computer 
Program 6 . 

The values of the eigenvalues obtained by the three 
methods are compared In Table I. In Table I are the eigen- 
values for the first differential equation, the difference 
equation with ten Intervals, together with the results of 
the computations, the number in parenthesis by an entry 
Indicates the number of Intervals used in the computation. 
The eigenvalues for forty Intervals are intermediate between 
those for ten and those for the differential equation. 



Differential 

Equation 


Difference 
Equation (10) 


Rayleigh 


Relaxation 


Dynamic 

Progranrolng 


4.44289 


4.42463 


4.4110 (25) 


4.42486(10) 


4.42442(10) 






4 . 44751(40) 


4.44611(40) 


4.43753(25) 


7.02482 


6.92714 


7.23029(40) 


6.92736(10) 


6 . 92717 ( 10 ) 


7.02482 


6.92714 


7.24707(40) 


6 . 92736 ( 10 ) 


6 . 92717 ( 10 ) 



Table I. 

Comparison of First Three Eigenvalues for Eq. (1.1) 



The three methods gave close agreement for the first 
eigenvalue but tended to differ on the next two. 

For the differential equation of higher order, only 
relaxation and dynamic programming were compared; these 
comparisons are shown in Table II using ten intervals. 
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Differential 

Equation 


Difference 
EquatlondO ) 


Relaxation 


Dynamic 

Programming 


19.739227 


19.577361 


19.60767 


19.57777 


49 . 3^80^10 


47.985220 


48.00555 


47.98473 


il9. 3^80^10 


47.985220 


48.02386 


47.98474 



Table II. 



Comparison of First Three Eigenvalues for Eq. (1.3) 

Computing times were similar, around twelve seconds 
for each. Computer Output 2 shows the results for the re- 
location method, and the program is Computer Program 2. 
Computer Output 9 shows the results for dynamic programming 
and the program is Computer Program 6. 



VII. DISCUSSION OP THE SOLUTION OF A 



PARTIAL DIFFERENTIAL EQUATION BY VARIOUS METHODS 

In Chapters III, IV, and V, three methods were developed 
for solving partial differential equations. These methods 
were used to obtain solutions of the following problems: 

V " - X - y) In A, (7.1) 

subject to the constraint 

U(x,y) =0 on 6A; (7.2) 

and 

V“U =8 in A, (7.3) 

subject to the constraints 

U(x,y) =0 on 6A 

— = 2(y -y"") for X = 0 (7.^<) 

9n^ 

^ = -2(y -y^) for x = 1 

3n^ 



= 2(x -x^) for y = 0 

3n^ 

— = -2(x -x^) for y = 0 

3n^ 

In this chapter the computation and numerical results are 
discussed . 
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In Galerkln’s method, an approximation for the function 
U(x,y) satisfying the constraint Eq. ( 7 > 2 ) was made by choos- 
ing suitable functions "^iCxjy), $2(x,y) and 'J>3(x,y) to satis- 
fy the constraint in Eq. ( 7 . 2 ). The approximation for U(x,y) 
was 

U(x,y) = Ci$i(x,y) + Cz^zU,y) + C3$3(x,y) ( 7 . 5 ) 

= (x -x^)(y -y^)(Ci +yxC2 +x^y^C3). 

The values of Ci, C2 and C3 were obtained by techniques 
described in Chapter V. Computer Output 3 shows the values 
of Cl, C2, C3 and the function U(x,y) at various points. 

In this method the region A was sub-dlvlded Into ten equal 
sub-intervals for each Independent variable for the inte- 
gration routine. A second approximation for U(x,y) was made 
for this method as 

U(x,y) = Ci$i + C2$2 + C3«I>3 ( 7 . 6 ) 

= (x -x^)(y -y^)(xCi + yCz + xyC3). 

The purpose of this was as follows. There generally is some 
skill and art involved in choosing the functions $^,..., 
well. In fact in the first choice is actually the desired 
function. The second set of functions was chosen so as to 
get some feel for the consequences of a poor choice of the 
The values of Cj, C^, C3 and the function at various 
points Is shown in Computer Output 10. 

The numerical solution of Eq. ( 7 . 1 ) by dynamic program- 
ming with the constraint of (7*2) yielded the values in 
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Computer Output 6. Again the region was divided as was 
done In Galerkln's method. 

In the Raylelgh-Ritz method, an approximation for the 
solution, U(x,y), satisfying Eq. (7.1) was made by choosing 
suitable functions $i(x,y), $ 2 (x,y) and '^’ 3 (x,y) to satisfy 
the constraint (7.2). The approximation for U(x,y) was 

U(x,y) = Ci$i + C 2^2 + C3^3 ( 7 . 7 ) 

= (x - x^)(y - y^)(Ci + xyC 2 + x^y^Ca). 

The values of Ci, C 2 , and C 3 were obtained by techniques 
described In Chapter IV. Computer Output 4 shows the val- 
ues of Cl, C 2 , C 3 and the values obtained for U(x,y) at 
various points. In this method It was found necessary to 
sub-divide the region into forty equal sub-intervals for 
each Independent variable in order to get satisfactory 
accuracy . 

The best approximation to U(x,y) was obtained by Galer- 
kln's method using Eq. (7.5), where $1 was the desired 
function. There was no error and the method was able to 
detect that this was the case. However, when Eq. (7.6) 
was used the maximum error was eight thousandths. The time 
required to solve the problem by this method was 0.57 
seconds . 

The problem was solved by dynamic programming in three 
seconds with accuracy to six digits. However, over half 
of the computer time was spent obtaining Inverses, which 
could have been fed as data from solving Eq. (6.1) In this 
particular case. 
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The Raylelgh-Rit z method required just under ten seconds 
and had a maximum error of two thousandths . 

Because of the time required and the accuracy obtained 
by Raylelgh-Rltz , only Galerkln's method and dynamic pro- 
gramming were used to solve Eq. (7*3)* 

Galerkln's method obtained the same values and accuracy 
In the solution of Eq. (7«3) as It did In the solution to 
Eq. (7.1) and took the same time. 

Dynamic programming obtained the same degree of accuracy 
In solving Eq. (7.3) as It did In solving Eq. (7.1) and 
took three seconds. The results are shown In Computer 
Output 8 and the program Is Computer Program 7 . 
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VIII. CONCLUSION 



The three methods considered for eigenvalue problems 
yielded satisfactory values for the eigenvalues and the 
eigenfunctions. Generally, the relaxation method seemed 
to be most satisfactory. It was straight-forward to pro- 
gram, and it was faster then Rayleigh-Ritz and dynamic 
programming. It converged rapidly even if step functions 
were used on several different tests figures, such as the 
L-shaped and triangular regions. 

The dynamic programming method converged in the same 
number of iterations as relaxation, but gave poorer esti- 
mates of the second and third eigenvalues. It of course 
was much more difficult to program and required more com- 
puting time due to the needed inverses. 

The Rayleigh-Ritz method seemed to have little to re- 
commend it due to the computer time required. It required 
much finer meshing in order to obtain a satisfactory accu- 
racy. The only advantage it had was that no iteration was 
required . 

Of the three methods considered in the solution of 
Poisson's and the biharmonlc equations Galerkln's method 
was the fastest and gave accuracy comparable to the Ray- 
leigh-Ritz method. It was also the simplest of the three 
methods used to program. 
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Dynamic programming gave the best accuracy generally, 
but It required more computer time than Galerkln's method. 

Raylelgh-Rltz had the same difficulties as it did in 
the eigenvalue problem and was considered of little use. 

While dynamic programming required more time for the 
solution of both eigenvalue problems and elliptic partial 
differential equations, it was very powerful. It only 
required a knowledge of the function U on the boundary. 

It can be extended to Irregular regions hh It obtained 
very good accuracy. Much of the time was spent computing 
the inverses. If the same points were used in solving 
several different problems, these inverses could be cal- 
culated once and thereafter entered as data, reducing the 
computer time greatly. 



50 



COMPUTER OUTPUT 1 . ( RE L AXATI ON ) 



EIGENFUNCTION=l 

PATH= 8 OMEGA= 4.424858 



X 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . ? 

0 . 9 

0 . 9 

0 . 9 



Y 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 



U(X, Y) 
0.09719044 
0.2508801 
0. 3094832 
0.2526796 
0.09705645 
0.2508801 
0.6512997 
0. 8061690 
0.6586339 
0. 2531579 
0.3094832 
0.8061690 
1.000000 
0.8169372 
0.3137754 
0.2526796 
0.6586339 
0. 8169372 
0.6658412 
0. 2551157 
0.09705645 
0.2531579 
0. 3137754 
0.2551157 
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COMPUTER OUTPUT 1 . ( REL AXAT ION ) 



EIGENFUNCTIGN=2 

PATH=11 OMEGA= 6. 927361 



X 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . I 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 

0 . 9 



Y 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 



U(X,Y ) 
0.1853049 
0.4686459 
0.6086538 
0.4958698 
0. 1895596 
0.3032387 
0. 8008109 
0. 9981067 
0.8116993 
0. 3099311 
0.001278793 
0.008393560 
0.01233941 
0 .008128799 
0.001685064 
-0.3083511 
-0. 8038315 
-0. 9953710 
-0.8089499 
-0.3103257 
-0.1921543 
-0.5CG6261 
-0. 6194299 
-0.5025631 
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CCMPUTER OUTPUT 1 . ( REL AXAT ION ) 



EIGENFUNCTION=3 

PATH=11 OMEGA= 6.927361 



X 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 

0 . 9 



Y 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 



U(X,Y) 
0. 1857997 
0.3045443 
0.002898388 
-0. 3070407 
-0. 1916554 
0.4894530 
0.8029368 
0.01102810 
-0. 8017028 
-0.4998162 
0.6086553 
0.9981196 
0.01235584 
-0. 9953650 
-0.6194320 
0.4950783 
0.8095817 
0.005508117 
-0. 8110784 
-0.5033802 
0. 1890559 
0.3086166 
5.9269C0'-05 
-0.3116446 
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COMPUTER OUTPUT 2 . ( REL AXATION ) 



EIGEN FUNCTION=l 

PATH=17 OMEGA= 19.60767 



X 



Y U(X,Y) 



0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 



0. 1011795 
0.257A610 
0.3156A26 
C.2600A89 
0.1008759 
0.2580318 
0. 6597930 
0.8128256 
0.6714125 
0.2608653 
0.3149233 
0.8084638 
1.000000 
0. 8271890 
0.3214785 
0. 2583916 
0.6639176 
0.8214287 
0.6778069 
0.2627081 
C. 1004035 
0.2579977 
0.3188558 
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CCMPUTER OUTPUT 2 . ( REL AX AT ION ) 



EIGENFUNCTI0N=2 

PATH=16 OMEGA= 48.00555 



X 



Y U(X,Y) 



0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 



0.1909801 
0.4904021 
0. 6059132 
0.4979762 
0.1925782 
C. 3037586 
0.7884125 
0.9826701 
0. 8096838 
0.3130895 
0.003429962 
0.01207077 
0.01897281 
0.01404283 
0.004551671 
-C. 3064085 
-0.7943833 
-0.9861545 
-0.8114642 
-0.3132687 
-0. 1937618 
-0.5022484 
-0.6227177 
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COMPUTER OUTPUT 2 . ( REL AX AT ION J 



EIGENFUNCTI0N=3 

PATH=19 OMEGA= 48.02386 



X 



Y U(X,Y) 



0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 



0.1851881 
0. 30 20115 
0.009329475 
-0.3049613 
-0. 1935259 
0.4791670 
0.7890626 
0.02798434 
-0.7933792 
-0.5036968 
0.5902148 
0.9785894 
0.03599443 
-0.9824019 
-0.6229793 
0.4810330 
0.7979080 
0.02484581 
-0.80 38256 
-0.5073704 
0.1851680 
0.3065956 
0.007537059 
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COMPUTER OUTPUT 3, (GAL ERKIN) 



VALUES OF C(I) ARE 



Cl = 1.000000 

C2 = 0 

C3 = 0 



X 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 



Y 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 



U(X,Y) 

0. 008099992 
0.01889999 
0.02249999 
0.01890001 

0.008100022 
0.01889999 
0.04409999 
0.05249999 
0.04410003 
0.01890005 
0.02249999 
0.05249999 
0. 06250000 
0.05250004 
0.02250007 
0.01890001 
0.04410003 
0.05250004 
0.04410006 
0. 01890007 
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COMPUTER OUTPUT 4. (RAYLEIGH) 



VALUES OF C(I) ARE 



Cl = 


1.120410 


C2 = 


-0. 5808935 


C3 = 


0.3500372 





X 




Y 


U(X,Y) 


0 


.25 


0 


.25 


0.03816108 


0 


.25 


0 


.50 


0.04937190 


0 


CM 

• 


0 


.75 


0.03599290 


0 


.50 


0 


.25 


0.04937190 


0 


.50 


0 


.50 


0.06231650 


0 


.50 


0 


.75 


0.04461559 


0 


.75 


0 


.25 


0.03599290 


0 


.75 


0 


.50 


0.04461559 


0 


.75 


0 


.75 


0.03179573 
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COMPUTER OUTPUT 5. (RAYLEIGH) 



EIGENFUNCTION=l 



THE 


VALUES 


OF C( I ) 


AR 


Cl 


= 


29.89044 




C2 


= 0. 


003808264 




C3 


= 


2.580980 






X 




Y 


0 


.25 


0 


• 


0 


.25 


0 


• 


0 


in 

• 


0 


• 


0 


.50 


0 


• 


0 


O 

LA 

• 


0 


• 


0 


.50 


0 


• 


0 


.75 


0 


• 


0 


.75 


0 


• 


0 


.75 


0 


• 



EGA= 4.447512 



U(X, Y) 

1.073552 
1.401157 
1.028184 
1. 431359 
1.868153 
1.370868 
1 .073485 
1.401070 
1.028116 



OM 

c 

25 

50 

75 

25 

50 

75 

25 

50 

75 
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COMPUTER OUTPUT 5. (RAYLEIGH) 



EIGENFUNCTION=2 


0M5GA= 


7.230292 


THE 


VALUES 


OF C(I) 


ARE 




Cl 


= 


9.664003 






C2 


= 


112.0784 






C3 


= 


-112.0294 








X 




Y 


U(X,Y) 


0 


.25 


0 


.25 


0.3401793 


0 


.25 


0 


.50 


1.766417 


0 


.25 


0 


.75 


2 .309447 


0 


.50 


0 


.25 


-0. 8598440 


0 


.50 


0 


. 50 


0.6040002 


0 


.50 


0 


.75 


1.765845 


0 


.75 


0 


.25 


-1.629948 


0 


.75 


0 


. 50 


-0.8604184 


0 


.75 


0 


.75 


0.3393202 
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COMPUTER OUTPUT 5. (RAYLEIGH) 



EIGENFUNCTION=3 OMEGA= 7.247070 



THE 


VALUES 


OF C( I ) , 


ARE 




Cl 


= 


9.656152 






C2 


= 


-112.4433 






C3 


= 


-111.6595 








X 




Y 


U( X.Y) 


0 


.25 


0 


.25 


-1.630179 


0 


.25 


0 


.50 


-0.8650630 


0 


.25 


0 


.75 


0.3325842 


0 


. 50 


0 


.25 


-0. 8558773 


0 


.50 


0 


.50 


0.6035087 


0 


.50 


0 


.75 


1.761141 


0 


.75 


0 


.25 


0.3463630 


0 


.75 


0 


.50 


1.770327 


0 


.75 


0 


.75 


2.309128 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

c 

0 

0 

0 

0 



COMPUTER OUTPUT 6. (DYNAMIC) 



X 



Y 



U(X,Y ) 



. 1 
• 1 
. 1 
. 1 
. 1 
. 3 
. 3 
. 3 
. 3 
. 3 
. 5 
. 5 
. 5 
• 5 
. 5 
. 7 
. 7 
. 7 
. 7 
. 7 
. S 
. 9 
. 9 
. 9 
. 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 



0.C08099854 
0.01889967 
0.02249958 
0. 01889964 
0.008099858 
0. 01889966 
0.04409914 
0.05249893 
0.04409913 
0.01889965 
0.02249958 
0.05249896 
0.06249879 
0. C5249900 
0.02249960 
C. C1889965 
0.04409916 
0.05249896 
0.04409920 
0.01889972 
0. C08099854 
0.01889966 
0.02249962 
0.01889972 
0.008099906 
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COMPUTER OUTPUT 7. (DYNAMIC) 



EIGEMFUNCTinN=l 



PATH= 


4 




OMEGA= 






4.423473 




X 






Y 




U( X,Y) 


0 


• 


1 


0 


• 


1 


0.09554338 


0 


• 


1 


0 


• 


3 


0.2500933 


0 


• 


1 


0 


• 


5 


C. 3091003 


0 


• 


1 


0 


• 


7 


0. 2500929 


0 


• 


1 


0 


• 


9 


0.09554338 


0 


• 


3 


0 


• 


1 


0.2500930 


0 


• 


3 


0 


• 


3 


0.6546414 


0 


• 


3 


0 


• 


5 


0.8090985 


0 


• 


3 


0 


• 


7 


0.6546412 


0 


• 


3 


0 


• 


9 


0.2500929 


0 


• 


5 


0 


• 


1 


C. 3091002 


0 


• 


5 


0 


• 


3 


C. 8090987 


0 


• 


5 


0 


• 


5 


1.000000 


0 


• 


5 


0 


• 


7 


0.8090994 


0 


• 


5 


0 


• 


9 


C. 3091004 


0 


• 


7 


0 


• 


1 


0. 2500930 


0 


• 


7 


0 


• 


3 


0.6546419 


0 


• 


7 


0 


• 


5 


0. 8090992 


0 


• 


7 


0 


• 


7 


0.6546427 


0 


• 


7 


0 


• 


9 


0.2500940 


0 


• 


S 


0 


• 


1 


0.09554315 


0 


• 


9 


0 


• 


3 


0. 2500930 


0 


• 


9 


0 


• 


5 


0.3091C10 
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COMPUTER OUTPUT 7. (DYNAMIC) 



EIGENFUNCTI0N=2 



PATH 


= 6 




OMEGA= 






6.927172 




X 






Y 




U (X,Y ) 


0 


• 


1 


0 


• 


1 


0. 1913275 


0 


• 


1 


0 


• 


3 


0.5004815 


0 


• 


1 


0 


• 


5 


0.6183150 


0 


• 


1 


0 


• 


7 


0.5004799 


0 


• 


1 


0 


• 


9 


0. 1913269 


0 


• 


3 


0 


• 


1 


0.3092604 


0 


• 


3 


0 


• 


3 


0.8089828 


0 


• 


3 


0 


• 


5 


0.9994567 


0 


• 


3 


0 


• 


7 


0.8089805 


0 


• 


3 


0 


• 


9 


0.3092589 


0 


• 


5 


0 


• 


1 


1 .959012'-06 


0 


• 


5 


0 


• 


3 


3.750642'-06 


0 


• 


5 


0 


• 


5 


2.407420'-06 


0 


• 


5 


0 


• 


7 


-4. 548319'-07 


0 


• 


5 


0 


• 


9 


-7.790408*-07 


0 


• 


7 


0 


• 


1 


- C. 3092566 


0 


• 


7 


0 


• 


3 


-0.8089775 


0 


• 


7 


0 


• 


5 


-0.9994535 


0 


• 


7 


0 


• 


7 


-0. 8089827 


0 


• 


7 


0 


• 


9 


-0.3092611 


0 


• 


9 


0 


• 


1 


-0. 1913257 


0 


• 


9 


0 


• 


3 


-0.5004784 


0 


• 


9 


0 


• 


5 


-0. 6183146 
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COMPUTER OUTPUT 7. (DYNAMIC) 



EIGENFU NCTI0N=3 

PATH=11 OMEGA= 6.927173 



X 



Y U(X,Y) 



0.1 
0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 



0.1913257 
C. 3092576 
-1.199204'-06 
-0. 3092595 
-0.1913273 
C. 5004785 
0.8089781 
-1. 8C1975'-06 
-0.8089815 
-0. 5004808 
0.6183135 
0.9994553 
1 .738089'-06 
-0.9994555 
-0.6183147 
0. 5004803 
C. 8089837 
4.765801'-06 
-0. 8089792 
-0. 5004812 
0.1913269 
0. 3092602 
2.622819'-06 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 



COMPUTER OUTPUT 8. (DYNAMIC) 



X 



Y 



U(X,Y) 



. 1 
. 1 
. 1 
. 1 
. 1 
. 3 
. 3 
. 3 
. 3 
. 3 
. 5 
. 5 
. 5 
. 5 
. 5 
. 7 

. 7 
. 7 
. 7 
. 7 
. 9 
. 9 
. 9 

. 9 
. 9 



0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 

0 . 1 

0 . 3 

0 . 5 

0 . 7 

0 . 9 



0.C08099731 
0.01889934 
0.02249917 
0.01889930 
0.008099727 
0.01889932 
0.04409828 
0.05249788 
0.04409826 
0.01889930 
0.02249917 
0. 05249792 
0.06249751 
0.05249795 
0.02249918 
0. 01889931 
0.04409831 
0.05249793 
0.04409834 
0. C1889938 
0 .008099716 
0. 01889932 
0.02249920 
0. C1889938 
0.008099772 
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COMPUTER OUTPUT 9. (DYNAMIC) 



EIGENFUNCTION=l 



PATH= 


3 




OMEGA= 






19.57639 




X 






Y 




U(X,Y) 


0 


• 


1 


0 


• 


1 


0.09549332 


0 


• 


1 


0 


• 


3 


0.2500034 


0 


• 


1 


0 


• 


5 


C. 3090197 


0 


• 


1 


0 


• 


7 


0.2500030 


0 


• 


1 


0 


• 


9 


0.09549332 


0 


• 


3 


0 


• 


1 


0.25C0032 


0 


• 


3 


0 


• 


3 


0.6545131 


0 


• 


3 


0 


• 


5 


0.8090189 


0 


• 


3 


0 


• 


7 


0.6545127 


0 


• 


3 


0 


• 


9 


0. 2500032 


0 


• 


5 


0 


• 


1 


0.3090197 


0 


• 


5 


0 


• 


3 


0.8090194 


0 


• 


5 


0 


• 


5 


1.000000 


0 


• 


5 


0 


• 


7 


C. 8090203 


0 


• 


5 


0 


• 


9 


0.3090200 


0 


• 


7 


0 


• 


1 


0. 2500033 


0 


• 


7 


0 


• 


3 


0.6545139 


0 


• 


7 


0 


• 


5 


0.8090203 


0 


• 


7 


0 


• 


7 


0.6545147 


0 


• 


7 


0 


• 


9 


0.2500044 


0 


• 


9 


0 


• 


1 


0.09549326 


0 


• 


9 


0 


• 


3 


0.2500034 


0 


• 


9 


0 


• 


5 


0.3090206 



67 



COMPUTER OUTPUT 9. (DYNAMIC) 
£IGENFUNCTION=2 



PATH= 


4 




OM5GA= 






47.98485 




X 






Y 




U( X, Y) 


0 


• 


1 


0 


• 


1 


C. 1911072 


0 


• 


1 


0 


• 


3 


0.5001713 


0 


• 


1 


0 


• 


5 


0.6181207 


0 


• 


1 


0 


• 


7 


C. 5001690 


0 


• 


1 


0 


• 


9 


0. 1911063 


0 


• 


3 


0 


• 


1 


C. 3091234 


0 


• 


3 


0 


• 


3 


0.8090394 


0 


• 


3 


0 


• 


5 


0.9998273 


0 


• 


3 


0 


• 


7 


0. 8090362 


0 


• 


3 


0 


• 


9 


0.3091207 


0 


• 


5 


0 


• 


1 


-9.384C42*-07 


0 


• 


5 


0 


• 


3 


8.968278*-06 


0 


• 


5 


0 


• 


5 


1.008165'-05 


0 


• 


c; 


0 


• 


7 


3. 054505*-06 


0 


• 


5 


0 


• 


9 


-4.620023'-06 


0 


• 


7 


0 


• 


1 


-C. 3091281 


0 


• 


7 


0 


• 


3 


-0.8090287 


0 


• 


7 


0 


• 


5 


-0.9998166 


0 


• 


7 


0 


• 


7 


-C. 8090367 


0 


• 


7 


0 


• 


9 


-0.3091323 


0 


• 


9 


.0 


• 


1 


-0. 1911175 


0 


• 


9 


0 


• 


3 


-0.5001757 


0 


• 


9 


0 


• 


5 


-0.6181255 
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COMPUTER OUTPUT 9. (DYNAMIC) 



EIGENFUNCTI0N=3 
PATH= 7 OME 

X 

0 . 1 
0 . 1 
0 . 1 
0 . 1 
0 . 1 
0 . 2 
0 . 3 

0 • 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 • 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 9 

0 . 9 

0 . 9 



47.98486 

Y U(X,Y) 

. 1 0.1911126 

. 3 0.3091245 

. 5 -9.687110'-07 

. 7 -0.3091258 

. 9 -0.1911120 

. 1 0.5001729 

. 3 0.8090308 

. 5 -2.217028'-06 

. 7 -0.8090366 

. 9 -0.5001741 

. 1 0.6181221 

. 3 0.9998204 

. 5 1.343255'-06 

. 7 -0.9998225 

. 9 -0.6181232 

. 1 0.5001751 

. 3 0.8090382 

. 5 5.371387'-06 

. 7 -0.8090329 

. 9 -0.5001737 

. 1 0.1911140 

. 3 0.3091279 

. 5 3.674680'-06 



GA 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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I 



1 




COMPUTER OUTPUT 1 0. ( GAL £RK I N ) 



VALUES OF C(I) ARE 



Cl = 1.596382 

C2 = 1.596392 

C3 = -2.598899 



X 



Y 



0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 1 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 3 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 5 

0 . 7 

0 . 7 

0 . 7 

0 . 7 

0 . 7 



0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

c . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 



1 

3 

5 

7 

9 

1 

3 

5 

7 

9 

1 

3 

5 

7 

9 

1 

3 

5 

7 

9 



U(X,Y) 

0.002375632 
0.01059511 
0.018627A9 
0.02069907 
0.01103618 
0.01059507 
0.03192532 
0.04658196 
0. 04633235 
0. 02294400 
0.01862741 
0. 04658184 
0.05916639 
0.05281764 
0. 02397245 
0.02069896 
0.04633217 
0.05281758 
0.04240138 
0. 01732974 
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COMPUTER PROGRAM 1. 



RELAXATION METHOD APPLIED TO 
V^U = -X^U, EQ. (1.1) 



B'GIN 

COMMENT FINDS f IGt NVALU;} S AND 21 GZi.'V; ^TC'.S BY 
P2LAXATI0N BAS’iD ^^l BAYLilGH’S FQPMULA ON AMY 
TYPS OF OHMS- IN. 

PORMAL PAR A MET' PS: 

GM HM’'GA 

0M2 OMZGA SOUAP'O 

DMC2 LAST VAl.DE OF OM-^GA SOUa'^^D 

ZNORM MAX. NORM 

Ml max. MUMPER OF X POSITIOMS 

'’I MAX. NUMBcR HP Y POSITIOMS 

HI MESH OF X VARIABLE 

H2 MESH OF Y VApIABLF 

M2 ARRAY OF LOW-R POSITIONS QP Y 

N3 ARRAY 0= U°P”R POSITIONS OF Y 

U,U1,Z UNKNOWN E I G^' MF ! )MC TI CMS 
PE POTENTIAL "N’RGY 

Ki KINETIC iNNRGY; 

R" 4 L KI.F:: ,0M2,0M0 2. Pa.ZS.OM.ZMORM; 

L H1,H2,C1 ,C2,C3,CAT, DOG, C AT L , DCG 1 ; 

INT"T,"R N,M,L,N1,P1; 

R’-^L H.H3,H4; 

0M2 :=5C00.0; 

N : = 1 ; 

£ :=0 .0C2 ; 

L : = 1 ; 

CATl:=50C0.0; 

R?fO(Nl,Pl,Hl,H2 ); 



H :=H1' m2; 
H3:=H1-^ -?; 
H4: =M- ■ ■-2 ; 

c- 1 M 

REAL ARRAY 
REAL ARRAY 



X( 0: :M1) 

Y ( 0 : : P 1 ) ; 



I NT 


g: r a 


Rr AY 


N2 , 


'^'3 ( 


J : :NT 


) ; 


REAL ARRAY Z,U 


,U1 


'( ^ i : 


: N 1 , 0 


: :P1) ; 


F'PR 


I : ='•• 


U NT T L 


Ml 


ro 






FOP 


J : =0 


I.INTI L 


PI 


00 








BEGIN 


U1 { I , 


J ) : 


=C . 








Z ( I , J ) 


• — r* 

• “ • V» 


f 










•-ND ; 












FOP 


!:=0 


!L 


N1 


DO 


X( I ) 


: = T "Hi; 


pop 


! : =0 


’JMTI L 


PI 


00 


Y( I ) 


:=I >H2; 


PGR 


I : =0 


UMTT L 


Ml 


on 


R’ AO ON (N2 (T 


POP 

Q • 


! : =0 


U \T I L 


Ml 


nn 

UV w 


R'-AD 


sOM(N3( I 


L ♦ 

IF 


L = 1 TH 














B“GIN 


L . • — ^ 












=0R ! : 


=r UNTTL 


Ml 


■IP 






FO® J : 


= 0 UNTIL 


= 1 


DO 






Z( I ,J) 


: = (X(I )- 


X(I 


) ' ^ 2 ) 


( Y( J )- 




r,n jn 


r • 

. 1 











IF 



iNO; 

L=2 T*J^ 
BEGIN 






=0R 

-OP 



;at 



Ml 

PI 



DU 

PQ 



I:=D UMTIL 
J:=n UNTIL 
BEGIN 

U( I,J):=ZC,J ); 

Z ( I , J ) : = ( X ( I ) -X ( I ) 

'MD; 

= K1 ;L : = 3 ;N : = 1 ;0M2 : =5000. 0 *, 



2 ) (Y ( J )-Y ( J ) V^2)t( 



K,' • 



5-X( ! ) ) ; 



71 



i 



I 



\ 




I 



I 



> 




■ 



I 



2)- (Y(J )-Y( J ) (0.5-Y( J) ) ; 



■11 TO D; 

'■fin ; 

IF I.. = 3 TH'M 

=OP I:=C '!\!T!L N1 
FOR j:='' M'lTIL PI DO 

B a T M 

IJKI , !) :=Z(I , J) ; 

Z( I,J): = (X(I)-X(I ) 

CAT r:' = K’’* ;i_ : = 4U':= l;OM2 : = 5000. 0; 

OT TO ?; 

CGWM-nf o:t and 

^ =K'*: =v. .i = ; 

PDF, I:=T UNTIL Nl-1 DO 
FOR J:=N2(I) UNTIL N3(I)-1 DO 
o"GIM 

PI:=F: + ( ( ( Z(I+; ,J)-Z( T ,J) )/H1 )^-2 ) = H 
+ (( (Z(I.J+1)-Z( T,J) )/H2)’'"2)?H; 

K-:=Kf:+(Z( !, J )' 2) H1^H2; 

* ND ; 

0M0 2:=‘jM2; 

0‘^2 ; =p^/K ■ ; DM: =SORT(D'^? ) ; 

TF (OM?>0MO2) (N>4C) TH'-’M 

F~OI^| 

WRf DT rnrJV 'p ) ; 

WRI TZ ( "N = " ) ; WF I TE ON ( N ) ; 

’JF,!T"‘(”C‘'''2 .f wr DN2 AP'. " ) ;WF!T?DN(0M02,0M2) ; 
on TO R; 

■2ND ; 

IP JRS ( 0M?-nv-i7 x'.'R TH-N 
BFPI N 

IPr,DNTFOL( ?) : 

WRiT„!(" IT” (•’ ");WPlT-(" •» ) ; 

WFIT'P(" COMPUT-'R OUTPUT 1 . ( R rL £ XAT 1 3N ) 

" ) ; 

WPTTCC ");WR!TP(" "); 
j j T n S T Z ' • =1 * 

wrtt":(»' . • X, -ighnfunction=».l-i) ; wRiTSc ••); 

INT = !’ LOST zr :=2; 

WRTT'-(« PATH=’SN,'’ 

WF I Ti: (" " ) ;WP ITT ( " ” ) ; 

WR I TF ( " X 

" U(x,Y)"); 

ViRiT* ( '' " ) ; 

= GR I:=l ST''P ? UNTIL M -1 00 
FOR J: = l ST ;p 2 UNTIL Pl-1 DO 
B?GIN IMT<=I iLOSIZ"': =2 ; 

WRIT"(" 'M DIV Nl, 

J DIV R-M P1,Z(I,J)); 

HPTTT(H II); 

•^ND ;~ND; 

IF (A3S(OM2-^’-T2)<-R.A ) ANF (L<4) TH'-N 
IF ABS('jM2-0'''02 X. RA THrN GO TO P ; 

N 1 ; 

COMN^NT Gl-T Z; 

FOR I:=l UHTTL Nl-1 DO 
BrGIN 

COP J:=N2(I)+1 UNTIL N3(I)-1 DO 

Z(! , J) : = (H4^ ( Z(I + 1 , J) + Z( I-l , J ) )+H3 •( Z(I t J+1 ) +Z( 1 1 J-1 ) ) ) 
/(2.‘-H3+2. HA-0M2 •( H>-'-2) ) ; 

"^ND; 

IF L>2 TM-'N 
B'GIN 

Cl :=OOG:=''OG1 :=0.C ; 

PQR I;=l "MTIL Nl-1 DO 
=0P J; = '‘'2(I) + 1 UNTIL N3(!)-l CO 
BOGIN 

C 1 : =C 1 + Z ( I t J ) N ; 

OOG:=OOG+Z(I,J) -U(I,J) H; 

OCGl :=D0G1 + Z( I , J) U1 ( I, J ) H; 



CN2GA = "t OM) ; 

Y " 



R=M Nl," 



GO TO e; 
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-MO ; 

C2: = COG/CAT ; C3 : =90G 1 /C;i T 1 ; 

^0“^ I:=l UNTIL Nl-1 O'l 

FOR J:=N2(!) + 1 UN Til N3(I)-1 DO 

Z (T , J ) ; = Z( I , J )-Cl-C2-' U( I , J )-C3-''Ul( I , J ) ; 

SND; 

D:ZN0F’M: = 0.0; 

FOR I:=0 until N1 DO 
>=CR J:=N2(!) UNTIL N3(!) O'T 
BEGIN 

ZS:=ABS(Z( I,J) ); 

1= ZS>ZMORM TH. M ZMORN: = ZS; 

END; 

FOR I:=0 UNTIL Ml DO 

FOP J:=N2(I) UNTIL N3 ( ! ) DO Z ( I , J ) : =Z ( I , J ) / ZNORM 
GO TO A; 

END; 

RiEND. 



73 



COMPUTER PROGRAM 2 



RELAXATION METHOD APPLIED TO 

V^U = A^U, EQ. (1.3) 



?■ GIN 

CrMM,;MT =!NDS ilG'NV&LUIS AND 5IG3NV~CTOPS BY 
P-LAXATIGM BASPD ON PAYL-IGH'S '=OR'^ULA FOR 
C! 4(ll)=W'*:2 U« 

= OP MAL PARAM-f T PS : 

OM OM'fPA 

GM2 OM3GA SQ'JARID 

CM02 LAST VALUP OF OMrGA SOUAR'^D 

ZNORM MAX. NORM 

H MISH SIZE OF X AND Y 

U,U1,Z UNKNOWN gIGENFUNCTIONS 

PE POTENTIAL ENERGY 

KS KINETIC ENERGY; 

REAL ARRAY Z , U ,U1 (C : : 1 4 , 0 : : 1 A ) ; 

REAL KE»P5.0M2.0M02,5RA, ZS.OM.ZNORM.H .CAT ,DOG , CATl ,DOGl; 
real Cl,C2tC3; 

INTEGER NtL; 

REAL ARRAY Y(Q: : 10); 

R'AL ARRAY X(0::10); 

LOGICAL SWITCH; 

REAL OMO; 

ON: =7000.0; 

SWITCH :=TRUE ; 

CAT:=CAT1:=2000.C; 

H:=0. l;0M2: =2000.0; N:=l; ERA: =0.003 ;L:=1; 

OM 2: =700000.0 ; 

^ P t =0 • C 0 5 * 

FOR 'l : =0 'UNTT L lo DO X ( I ) : =Y ( I ) : = I<=H ; 

‘=0R I:=0 UNTIL 14 DO 

>=0R J:=C UNTIL 14 DO Z ( I , J ) : =U ( I , J ) : =U 1 ( I . J ) : =0 . 0 ; 

BEGIN 

PROCEDURE NORMA; 

begin ZN0RM:=0.0; 

FOR I; =3 until H DO 
FOR J;=3 UNTIL 11 00 

BEGIN ZS:=ABS( Z< I , J ) ) ; 

1= ZS>ZNO.RM THLN ZNOR‘^: = ZS;r;ND; 

FOR I:=0 UNTIL 14 DO 

FOR J:=0 UNTIL 14 DO Z ( I , J ) : =Z ( T , J ) / ZNORM; 

GO TO A; 

END NORMA; 

PPCCEDURg FIJNCT; 

BEGIN 

IF L = 1 TH'EfJ 
BEGIN L:=2; 

=0R I:=0 UNTIL 10 DO 

FOR J:=0 UNTIL 10 DO 

Z( 1+2, J + 2) :=( Xn ) - X(! )*>««2)?^ ( Y( J ) -YU )*t2 ) ; 

NORMA ; END; 

IF L=2 TH«N 

BEGIN l.: = 3 ;C AT :=<E ;N: = l;0M2 : = 2000.0; 

CM: =7000.0 ; 

FOR I:=0 until 14 DO FOR J:=0 UNTIL 14 DO U (I , J ) : =Z ( I , J ) ; 

FOR I:=0 UNTIL 14 DO FOR J : =0 UNTIL 14 DO Z(I,J):=0.0; 

FOR I:=0 UNTIL 13 DO 

FOR J:=0 UNTIL 10 DO 

Z(I+2,J + 2): = (X(I)-X(I )*«2)«( YU )-Y( J ) ^i*2) M . 5-X( I ) ) ; 
NORMA; END; 

IF L=3 THEN 

BEGIN L: =^;CATl:=XE ;N : = 1 ; OM 2: = 2000 . 0 ; 

OM: =7000.0; 

FOP l:=0 UNTIL 14 DO =0P. J : =0 UNTIL 14 00 U1 ( I , J ) : = Z ( I ,J ) 
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<=0» I:=n UNTIL 14 DO FOR J:=0 UNTIL 14 DO Z(l7J):=0.0 
I ; = n UMTI L 10 DO 
FOR J:=0 UNTIL 10 DO 

Z(I+2 ,J + 2):=(X( I )-xn )a-*2)=«'(Y(J)-Y( J)**2)-^( .5-Y(J ) ) 
NORMA;FND; 

2ND PUnCT; 

IP SWITCH ^H'fN BrGIM S WI T CH r = FALS 5 ; FUNCT ; END *, 

COMMENT GET FE AND KE ; 

A: PE:=KE:=0.r ; 

FQp I: =2 UNT^L 12 DO 
FOR J:=2 UNTIL 12 DO 

BEGIN PJ:=?E+(Z( I+1,J ) -4. *Z (I , J ) +Z ( I- 1 , J ) +Z ( I , J+ 1 ) 

+ Z(I » J-1) )T.^2/(H-i=H); 

•<E:=KE + (Z( I,J *«2;EMD; 

OMO:=OM; 

0M2:=?c/KR ;OM:=SQRT( 0M2) ; 

IF (CM>OMO) OP. (N>60) THEN 

BEGIN WRTT.'C'NOT CTNVERGING 
WRITEC ") tWRITCf ''MtZNQRM,OMO,OM'’) ; 

WRITS (M,ZMURM,OMOt DM) ;WRITE (" " ) ; 

FOR I:=l UNTIL II DO 

WRITc(Z( 3t I ) ,Z( 5. I ) ,Z(7, 1) , Z(9» I ) »Z (11 ) ); 

GO TO R:”ND; 

IF ABS(CMO-OMXFRA THEN 
BEGIN lOCONTROL ( 3 ) ; 

WRITEC WRITE (" ”); WRITE ('* "); 

WRITSC COMPUTER OUTPUT 2.”» 

•• (R"LAXATI0M) •’ ) ; 

WRITEC ") ;WRIT£( " '• ) ; 

INTFnLDSIZE:=l; 

WRITEC’ EIGENFUNCTI0N=»,L-1) ; WRITE( ” ") 

INTFISLDSI Z"'’:=2; 

WRITEC PATH = 'SN," OMEGA = ".CM); 

WRITEC " ) ;WRITEC '' ) *, 

WRITEC X Y 'S 

” U( X. YCM ;WRI TS( ” •’); WRITEC* ") ; 

POR I:=l STEP 2 UN'^R 9 DO 
FOR J;=l STEP 2 UNTIL 9 DO 
BEGIN INTFiaOSIZ5;=2; 

WRITEC ",I DIV 10, ".'M R£M 10, 

" ",J DIV 10,*'.", J REM 10 , Z( 1+2 , J+2 ) ) ; 

WRITS (" " ) ; 

END; 

IF L<4 THEN FUNCT ELSE GO TO R; 

END; 

N: =N+1 ; 

COMMENT GET M'^W Z; 

FOR I:=3 UNTIL 11 DO 
PQR J:=3 UNTIL 1 1 DO 

Z( I , J) : = ( B.-.-Z( I+l , J ) + 8CZ (I -1 , J) + 8X- Z( I , J+1 ) 

+8.-Z( I , J-1 )-Z( !+2,J)-Z(I,J-2)-2.*Z(I+l,J+l) 
-2.1=Z(I-1 ,J-1 )-2.';-Z{I-l ,J+1 )-2.--<3Z{I + l,J-l ) 

-Z( 1-2, J )-Z (I , J + 2) ) /( 20.-0M2^?H“Hr H*-H) ; 

FOR I:=2 UNTIL 12 DO 

BEGIN Z(0,I ) ;=-Z(4,I) ; Z(14,I) :=-Z(10,I ); 

Z( 1, 1 ) ;=-Z (3, I ) ;Z (13, 1 > :=-Z( 11 ,! ) ; 

Z(I ,0):=-Z(!,4);Z(I,14): = -Z(I,lC ); 
Z(!,1):=-Z(I,3);Z(I,13) :=-Z(I ,1 1) ;cND; 

IF L>2 THEN 

BEGIN Cl : =00G:=D0G1 :=0 .0 ; 

FOR I: =3 UNTIL 11 DO 

FOR J;=3 UNTIL 11 DO 

BEGIN Cl :=C1 + Z( I , J) *H-^H; 

DOG:=DOG+Z( I , J )-''U( I , J )=i^H'^H; 

DOGl:=DOGl + Z( I , J )’^U1 ( I , J ) -H^^H ; EN D ; 

C2:=D0G/CAT;C3 :=DCG1/CAT1 ; 

FOR I:=3 UNTIL 11 DO 

FOR J:=3 UNTIL 11 DO 

Z( I , J) :=Z( I , J)-C 1-C2'^U(I , J) -C3I-U1 (I , J) ; 

FOR I: = 2 LV'TIL 12 DO 

BEGIN 

Z(0,1 ) :=-Z(4,I ) ;Z ( 14, 1 ) :=-Z(10.I ) ; Z ( 1 , 1 ) : =-Z ( 3 , 1 ) 
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ilJOL 



Z(13,I):=-Z(11,!);Z(I,0):=-Z(I,4);Z(It 14) i=-Z( I , 10) 
Z(I ,1) :=-Z(I,3);Z(I,13):=-Z(I,ll);=ND; 

END; 

IF L<5 THEN NORMA; 

ND; 

;FND. 



COMPUTER PROGRAM 3 



GALERKIN'S METHOD APPLIED TO 
V^U = 2(x^+y^-x-y) , EQ. (1.5) 



BEGIN 

CCMMcNT SOLVES PARTIAL DIFFERENTIAL EQUATIONS QP 
THE FORM L(U)=F BY GALERKIN'S METHOD ON ANY DOMAIN. 
FORMAL PARAMETERS: 

N NUMBER OF EQUATIONS 

M2 MAX. NUMBER OF X OOSITICNS 

M3 MAX. NUMBER OF Y POSITIONS 

N2 ARRAY OF LOWER Y POSITIONS 

N3 ARRAY 0= UPPER Y POSITIONS 

F KNOWN FUNCTION 

U UNKNOWN FUNCTION 

C(I) ESTIMATING FUNCTIONS 
OKI) L(SSTIMATING FUNCTIONS) 

A(I) COcF. OF ESTIMATING FUNCTIONS; 

REAL H,E,DOG,CAT,DEV,AMUL,T, OOM; 

REAL Hi; 

REAL H2.H3; 

INTEGER M2, M3; 

INTEGER N,M,R ; 

E : =0 .COOOO 1 ; R ; =0 ; DEV : = 1 . 0 ; 

READ (N, M2, M3 ,H2,H3); 

Hl:=H2-«‘H3; 

BEGIN 

REAL ARRAY X(0::M2); 

REAL ARRAY Y(0: :M3) ; 

REAL ARRAY S( 1 : : N, 1: : N + 1 ) ; 

REAL ARRAY A,C(l: :M) ; 

REAL ARRAY F . U ( 0 : : M2 , 0 : : M 3 ) ; 

REAL ARRAY 0 , 01 ( 1 : :N ,0 : : M2 ,0 : : M3 ) ; 

INTEGER ARRAY M2 , N3( 0 : : M2 ) ; 

FCR l:=0 UNTIL M2 DO X(I): = I’^H2; 

FOR I:=0 UNTIL M3 00 Y(I): = I?=H3; 

FOP I:=0 UNTIL M2 DO RcADCN ( N2 ( I ) ) ; 

FCR I:=0 UNTIL M2 DO RE ADON ( M3 ( I ) ) ; 

FOP I :=0 UNTI L M2 DO 

FCR j:=0 UNTIL M3 DO 
BEGIN 



F(I ,J) :=0. O; 
U(I ,J ):=0.0; 
FOR L:=1 until 



N 



D(L,I , J) :=01(L, I , J) :=0.0; 

END; 

BEGIN COMMENT SOLVE HERE; 

PRCCEOURE WRITER; 

BEGIN WRITSC'S INGULAR" ) ;G0 TO Q; 
END WRIT':®; 

PROCEDURE PE‘5F0R.M; 

BEGIN COMMENT GET U(I,J) AND A( I) ; 



") ; WRITE (" " ) 



WRITEC ");WRIT£(" ");WRIT£(" 
WRITE( •' 

WRITE ( " 

WRI TE ( " 

WRITEC* "); 

FCR I:=l UNTIL N DO 
BEGIN 

Ad) : = B( I, N+1 ) ; 
INTF!5LDSIZE:=l; 

WR I TE ( ” 

END; 

FOR J:=0 UNTIL M3-1 DO 
BEGIN 

FCR I; =0 UNTIL M2-1 CO 
BEGIN DOG: =0.0; 



") ; 

COMPUTER OUTPUT 3.(GALERKIN) 
VALUES OF C(I) ARE") ;WRITE(" "); 



C", I ," = ",A(I) ); WRITEC* ••); 



'*) 
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■^CiP L: = 1 until N DO DOG : =DOG+ A ( L ) + D ( L , I ♦ J ) ; 

U( I ,J):=DJG; 
f n'D;FND; 

WRITEC’ ");WRIT~(» ••);WRITE(” ; 

WRITEC X Y ", 

"U( X,Y) ") ; 

WRIT^= (" " ) ;WRITE(" " ) ; 

FOR l:=l ST*P 2 UNTIL M2-1 DO 
FCF j:=l STEP 2 UNTIL M3-1 DO 
BEGIN INTFIELDSIZE:=2 ; 

WRITcC ",I DIV M2,".'M REM M2, 

" ",J DIV REM M3,U(I,J)); 

WRITE!" "); 

END ; 

GC TO Q; 

END PER'^ORM; 

PROCEDURE SWITCH; 

BEGIN CCMMZ-MT SOLVES FOR A(I) COcFF; 

INTEGER Ml,vn_; 

Ml;=N;Ml:=N+l; 

FOR l:=l UNTIL Ml DO B ( I . Ml + 1 ) :=C ( I ) ; 

FOP K;=1 UNTIL Ml DO 

BEGIN IF K=Ml them GO TO P; 

POR I:=K+1 UNTIL ^1 DO 

BEGIN IF ABS( B( K, K) XABS ( B( I ,K) ) THEN 
BEGIN R:=R+l; 

FOR J:=l UNTIL N1 DO 

BEGIN T: = 3(I,J);B(I,J):=B(K,J);B(K,J): = T ;END 

5ND;cNO ; 

P:IF ASS(B(K,K))<F THEN WRITER ELSE 

BEGI N D£V:=DEV> B( K ,K) ; DGM: = B(K, K ) ; 

FOP J; = l UNTIL N1 DO B ( K , J ) : = B ( K , J ) /DOM; 

FOR l:=l UNTIL Ml CO 

BEGIN AMUL:=B( I ,K); 

IF I=K THEM ELSE 
BEGIN 

FOR J:=l UNTIL N1 DO 

BEGIN B( I ,J ) ;=B( I , JI-AMULI'B! K.J) ; 
£ND;END; END;ENO;£ND; 

PERPOF m; 

END SWITCH; 

COMMENT GUESSED pUNCTICNS D(I) HERE; 

FOR J:=0 UNTIL M2 DO 
FOR K;=N2(J) UNTIL N3(J) DO 
BEGIN 

F( J,K) :=2.=^(X( J)>«''^2 + Y(K):«-4:2-Y(K)-X( J) ); 

D( 1 tJ.K ); = ( X( J )-X( J )’ *2 )^( Y( K )- Y( K X*2) ; 

C(2 ,J,K ): = (X( J)*^2-X( J )^«3 ) ’'• ( Y ( K ) ^ 2-Y ( K ) 2 ) ; 

Cl(l,J,K);=2.*(X(J)»t2 + Y(K)>^?2-X(J)-Y(K)); 

Cl (2. J.K ) : = (2.-6.^X( J ) )M Y( K 2-Y(K ) Y.--^3) 

+ (2.-6,-P-Y(K) )- (X ( J)^=P2-X ( J )’^^5) ; 

END; 

COMMENT SOLVE FOR INTEGRALS F^D AND LD^D; 

FOR I : =1 U^!TI L N DO 
BEGIN D0G:=0.0; 

FOR J:=l UNTIL M2 DO 

FOP K:=N2(J)+1 UNTIL N3 ( J ) DO 

C0G: = D0G+( F( j ,K )=» D( I, j ,K ) t'Hl ) ; 

Cd ) :=DOG; 

FOP J:=l UNTIL N DO 
BEGIN CAT: =0.0; 

FOR K:=l UNTIL M2 DO 

FOP L:=N2(K)+1 until N3(K) DO 

CAT: = CAT + ( Cl( J,K, L )«D( I,K,L )’^H1) ; 

B(I ,J) :=CAT; 

END ; 

IF I=N THEN SWITCH; 

END; 

END ;SND; 

0:END. 
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COMPUTER PROGRAM . 
RAYLEIGH-RITZ METHOD APPLIED TO 
V^U = 2(x^+y^-x-y) , EQ. (1.1) 



8' GIN' 

COMVi.NT SOLVES PfRTUL D !'='='= R’: NT! AL ‘EQUATIONS BY 
RAYL^IGH-PITZ M' TMDD. 

FORV' L '^ARAN3T’RS : 

H2 OX S^.^CING 

H3 OY SPACING 

N NUMBER .3Q'J/TIONS 

F KNOWN FUNCTION 

M2 UPPTR LIMIT OF X 

V3 UFP".R LIMIT CP Y 

u unknown function 

N3 UPP'R BOUNDS ON Y 

N2 LOW”P BOUNDS ON Y 

A(I) CO'^ = . QP PSTIMATIMP, FUNCTIONS 

D(I) pSTTMATING PU^ICTTOMS; 

R5AL H,i,DPG,CAY,cr V,AMIJL,T,0GM,H1 ,H2 ,H3; 

INT"G''P 

p : =C . 00000 1 ; F : =0 ; r IV ; = 1 . u ; 

RPAD(.N,M2,M3,H2,H3); 



Hi: =H2 -'W3 ; 

BCG IN 

RFAL ARP-JY X(C::M2); 

RSAL t. RRmY Y (P: : M3) ; 

REAL ARR5\Y G ( 1 : ; N , 1 : : N+1 ) ; 

REAL ARRAY A,C ( 1: ;N ) ; 

RPAL ARp./Y =,U(0::M2,C::M3) ; 

RTAL APPAY D ,D 1 ,D2 ( 1 : : N , p : : M2 ,C : : M3 ) ; 

INTEGER aRpAY M2,N3(p : :M2 ) ; 

FCR I : =0 UNT I L M2 X ( I ) : = I ^"^2 ; 

FOP I:=0 UNTIL M3 DO Y(I):=I'H3; 

FOP I:=0 UNTIL M2 DO F.f A DOM ( N 2 ( I ) ) ; 

cnq I ; =M M2 DO P "’aDP'NI N3 ( I ) ) ; 

POP I:=0 UNTIL M2 DO 

FOR J:=0 until M3 DO 

P.:-GIN 

F( I ,J) :=iJ( I ,J ) : =0. 0; 

=0P L:=l UNTIL N DO 

D(L , J) :=D1 ( L ,I , J) :=D2( L, I , J ) : =0 .0; 

TMD I 

FOR I:=l UNTIL N 00 

FOR J:=l UNTIL N+1 DO B(I,J):=C.P; 

B? GTN 

C.DMMPMT SOIV. H'lR-'; 

RROCrTUPf WR! T. R ; 

BFGtm WR it "S INGULAR" ) ;go to 0; 

INC WRIT-R; 

PPOC'DIJR- F'lPPQRM; 

^ ^ G 1 M 

CCmm;nt git U(I,J) AMD Ad); 

WRJTFI" ");writ:-(" ");wr!t:c* 'M; 

WRTT“(" COMPUTER OUTPUT 4. (PAYLf 

WR! T’ ( " " ) ; write ( »' »» ) ; WR I T I ( " " ) ; 

WPIT"(" VALUES OF C(I) AR i" ) ; WR I TE ( •' " 

=0R I:=l UNTIL N CO 
6;GI N 

A ( T ) :=P( I ,M+1) ; 

INTFI”LDSIZl:=l; 
wRj t;- ( '• 

' MC ; 

FOR I:=l UNTTL 'd- 1 
= PR j:=i until N3-1 
G"GIN DOG: = 0.0; 

for L:=1 until N DO DOG: = DOG+A ( L )-0 ( L , I t J ) ; 



C, = A( I ) ) ;WRIT”( " " ) ; 



DO 

CD 



IGH ) " ) 
); 
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U(! ,J) :=D">G; 

■-ND ; 

WRIT“(" 'M'.WRITrC' '• ) ;V-^RIT.:('' "); 

wR! n ( " X Y , 

'• u ( X ♦ Y ) " ) ; 
wRlTrC ") ;VJDIT5 C " ) ; 

FOP. I:=0 ST“t^ 10 UNTTL M2 DO 
POP J:=J STrP 10 UNTIL M3 00 
BEGIN INT'^IRLDSIZ.-:=2; 

WRITEC "t! Div M 2 ,''.'M-nr;r R"M M 2 , 

•• ",J riV M2, Jr 100 PEi^ M3,U(I,J)); 

WRITP(" '•); 

“NC ; 

GO TO 0; 

END PFPFORM; 

PROCECUP.P SWITCH; 

BEGIN 

comment SHLV: S fop A(D CC^F = ; 

INT-G'R 

Ml ; =N; t\l ; = N+1 ; 

FOP I:=l UNTIL Ml DO S ( I , Ml + 1 ) : =C ( I ) ; 

FOP K:=1 until Ml DO 
BEGIN 

IP K=M1 THFM. GC TO P; 

FOF I:=K+1 UNTIL Ml 00 
BEGIN 

IF ABS( B( K,K) )<ABS(B( I ,K) ) THEN 
BEGIN 
R: =? + l ; 

FOP J;=l UNTIL N1 DO 

BEGIN T; = B( I, J ) ;B( I ,J ) ;=B(K,J) ;B( K,J) :=T;3ND; 

~ N 0 * ’ * r * n * 

P:IF ABS( Bfk ;k ) )<r THEN WF.IT^p f- L SE 
B'GIN 

0~V:=O^V«B(K,K) ;rOM:=P(K,X) ; 

=0R J: = l until N1 00 B ( K , J ) : =B ( K , J ) ZOOM ; 

PQR I:=l UNTIL ’■'1 DO 
B”GIN 

A‘-'iJL :=F( I ,K ) ; 

IF I=K THrN 'LS^ 

BT’OIN 

FGP J:=l UNTIL N1 00 
BEGIN 

B(I ,J) :=B(I ,J) -AMUL' PIK, J) ; 

END ; END; 7N0 ; HND;pnd; 

PEPPORM; 

'^ND SWITCH: 

COMMENT GU SSiO FUNCTIONS 0(1) H7 ?p ; 

POP j:=0 UNTIL M2 DO 
FOR K:=N2(J) UNTIL N3(J) 00 
BEGIN 

F ( J ,K) :=2. M X( J)^ 2+Y -X ( J )-Y ( K ) ) ; 

D( 1, J ,K ) : = ( X( J )-X (J ) ' ‘.E) M Y( K) -Y( X ) v=r. 2) ; 

D(2,J,K):=X(J) Y(X) 

•^OtUtKIr^XU) YC<) 0 ( 2 ,J,K); 

D1(1,J,K ) ;=(l.-2. X(J)) ( Y(K ) -Y(K )^^T2) ; 

01 ( 2, J,K ) ; = (2. X( J) -3, -X ( J )-. -'2 ) ( Y ( K ) •-' 2- Y ( K ) ■'" 3 ) ; 

D1(3,J,K) :=(3. X( J) >•.<2-4. -X(J) -3)- (Y(t<)>t3-Y(K)i'»4); 

1^2 ( 1,J ,K ): = (1.-2. Y (K ) ) ( X ( J ) - X ( J ) ■ 2 ) ; 

D2( 2, J ,K) : = (2. V( K) -3.= Y(K ) ^'-2 ) ( X ( J ) " ' 2-X ( J ) >- 3 ) ; 

D 2 ( 3 , J,K ) : = ( 3 . Y(K ‘ 2 - 4 . ty( K)-' 3 ) ( X ( J ) • - 3 -X ( J ) - 4 ) ; 

’ND; 

SOLVE =nR IM^^.GRa.LS NOW; 

FOP I:=l IPiTIL N DC 
BEGIN 
DOG:=C . C; 

FOR J:=l UNTIL M2 DO 

K:=M2(J)+1 until N3(J) DO 
DOF : =DOG- ( = ( J , K ) * D ( I , J , K ) >'H1 ) ; 

C( T ) : = DCG; 

-NO; 

FOP !:=I UNTIL N CO 
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POR J:=l 

QV T M 

- 0 " ‘k : 

!=0P L: 

en , j) 

HI ; 
'-' 10 ; 
swi TO’; 
jN 0 ;''ND; 
OtiMO. 



nr 'TO N 00 

= 1 'HTIL M2 DO 
= r'2( *< ) +1 UNTIL N3( K) DO 

: = 3(I, J ) + (Dl ( I, K,L )=?:C1( J,K, L )+02( I,K,L ) 02( J,K,L) 
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COMPUTER PROGRAM 5. 
RAYLEIGH-RITZ METHOD APPLIED TO 
V^U = -X^U, EQ. (1.1) 



BEGIN 

COMMENT SOLVES EIGENVALUE PROBLEMS BY RAYLEIGH-RITZ 
METHOD USING TRE02 AND TGL2 FROM NUMERISCHE 
MATHEMATIK 11, PP. 181-195 AND PP. 293-306(1968) BY MARTIN 
REINSCHjBOWDLER ,HILARY, AND WILKINSON. 

FORMAL PARAMETERS: 

N NUMBER OF EIGENVALUES 

AND THE NUMBER OF EQUATIONS USED 
HX MESH OF X VARIABLE 

HY MESH OF Y VARIABLE 

NPX NUMBER OF X POSITIONS 

NPY NUMBER OF Y POSITIONS 

U(I) UNKNOWN EIGENFUNCTIONS 
C(I) ESTIMATING FUNCTIONS USED 
R(I) COEF. OF ESTIMATING FUNCTIONS 
02(1) EIGENVALUES OF MATRIX A 
DKI) EIGENVALUES OF MATRIX B 
ZKI) EIGENVECTORS OF MATRIX B 
Z2(I) EIGENVECTORS OF MATRIX E=T'*A*T 
D(I) USED IN TRANSFORMATION C=T'=f=D; 

INTEGER NX, NY, N, NPX, NPY; 

REAL HX,HY,EPS,DOG; 

LOGICAL T; 

PROCEDURE TRED2( INTEGER VALUE N;REAL ARRAY D,E(*); 

REAL ARRAY A , Z ( * , ’f' ) : R E AL VALUE PU ) : 

COMMENT REDUCES SYMMETRIC MATRIX TO TRIOIG. FORM USING 
HOUSEHOLDER, S REDUCTION. 

FORMAL PARAMETERS; 

N ORDER OF SYMMETRIC MATRIX A 
0 DIAGONAL OF RESULTS 
E SUB-DIAGQNAL OF RESULTS; 

BEGIN 

INTEGER I,J,K,L; 

REAL F,G,H,HH,TOL ; 

tol;=pu/epsilon; 

FOR l;=l UNTIL N DO 

FOR J: = l UNTIL I DO Z ( I , J ) ; =A { I , J ) : 

FOR Il;=N STEP -1 UNTIL 2 DO 
BEGIN I:=Il; 

L;=I-2;F;=Z(I ,1-1) ;G:=0.0; 

FOR K;=l UNTIL L DO G : =G+ Z ( I , K ) «Z ( I , K ) ; 

H:=G+F*F; 

IF G <= TOL THEN 

BEGIN E( I ) :=F ;H:=0.0;GO TO SKIP;ENO; 
l;=L+1 ; 

G:=E(I);=IF F>= 0.0 THEN -SQRT(H) ELSE SQRT(H); 
h:=h-f*G; Z( 1,1-1) ;=F-G;F;=0.0; 

FOR J;=l UNTIL L DO 
BEGIN 

Z( J,I ) ;=Z(I ,J)/H;G:=0.0; 

FOR k;=i until j do G ; =G+Z ( J , K ) *Z ( I , K ) ; 

FOR K;=J+1 UNTIL L DO G: =G+Z ( K, J ) *Z ( I , K ) ; 

E(J) ;=G/h;F;=f+Gyz(J,I) ; 

END; 

HH:=F/(H+H) ; 

FOR J;=l UNTI L L DO 
BEGIN 

F;=Z( I , J) ;G: = E( J) ; = E( J)-hh=!^f; 

FOR K;=1 until j do 

Z( J,K) :=Z (J ,K)-F*E(K)-6«Z( I,K) ; 

END; 

SKIP:d(I ) ;=H; 
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END; 

D(l) :=E(1) ;=0.0; 

FOR I :=1 UNTIL N DO 
BEGIN L:=I-l; 

IF 0(1) -.= 0.0 THEN 
FOR J;=l UNTIL L DO 
BEGIN G;=0.0; 

FOR K:=1 until L DO G : =G+Z ( I , K ) *Z ( K, J ) ; 

FOR K:=l UNTIL L DO Z ( K , J ) : =Z ( K, J ) -G* Z (K , I ) 
END; 

D(I ) ; = Z(I , I ) ; Z(I , I ) : = l.O; 

FOR J;=l UNTIL L DO Z ( I , J ) : =Z ( J , I ) : =0.0 ; 

END; 

END TRED2: 

PROCEDURE TQL2(INTEGER VALUE N;REAL ARRAY D,E(*); 

REAL ARRAY Z ( >>= , =5^ ) ; L OG I C AL RESULT FAIL): 

COMMENT FINDS EIGENVALUES AND EIGENVECTORS. 

FORMAL PARAMETERS: 

N ORDER TRIDIAGONAL MATRIX 
D DIAGONAL ELEMENTS 
E SUB-DIAGONAL ELEMENTS 

Z MATRIX OF HOUSEHOLDER TRANSFORMATION; 

BEGIN 

INTEGER I,J,K,L,M; 

REAL B,C,F ,G,H,P, R, S; 

fail;=false; 

FOR l:=2 UNTIL N DO E(I-1):=E(I); 

E(N) :=B:=F:=0.0; 

FOR Ll:=l UNTIL N DO 
BEGIN L:=L1:J:=0; 

H:=EPSILON=i=(ABS(D(L) ) +ABS (E(L) ) ) ; 

IF B<H THEN B :=H; 

COMMENT LOOK FOR SMALL SUB-DIAGONAL ELEMENT; 

FOR Ml :=L UNTIL N DO 
BEGIN M:=Ml; 

IF ABS(E(M)X=R THEN GO TO CONT;END; 

CONT-.IF M=L THEN GO TO ROOT; 

NEXTITriF J=30 THEN 

BEGIN FAIL:=TRUE;G0 TO FIN; END; 

J:=J + 1 ; 

P;=(D( L+1 )-D( L) ) / (2.*E(L) ) ; 

R: = SQRT( P**2+ 1.0 ) ; 

H;=D(L)-E(L)/{ IF P<0.0 THEN P-R ELSE P+R); 

FOR I:=L UNTIL N DO D(I):=D(I)-H; 

F;=F+H; 

COMMENT QL TRANSFORMATION; 

P;=D(M) ;C:=1.0:S:=0.0; 

FOR l:=M-l STEP -1 UNTIL L DO 
BEGIN 

G:=C*E( I) ;h:=C*P; 

IF ABS(P) >= ABS(E(I ) ) THEN 
BEGIN 

C;=E( I ) / P;R; = SQRT ( C^=«' 2+1.0 ) ; 

E( I+l ) : = S=f=P*R;S:=C/R;C;=1.0/R; 

END ELSE 
BEGIN 

C: = P/E( I ) ;R:=SQRT (C**2+1.0) ; 

E(i+i ) : = s=i‘E(l )*r;s;=1.0/r;C;=C/r; 
end; 

P:=C*D(I)-S*G; 

D( I + l ) :=H + S=«'( C'?'G + S*D(I ) ) ; 

COMMENT FORM VECTOR; 

FOR K;=l UNTIL N DO 
BEGIN 

H;=Z( K, I+l) :Z(K,I+1 );=S*Z(K, I )+C*H; 
Z(K,I ) :=C«Z(K, I )-S*H; 

END;END; 

E(L) : = S*P;0(L ) :=C’f=P; 

IF ABS(E(L)) > B THEN GO TO NEXTIT; 

ROOT :D(L) :=D(L)+F; 

END; 

COMMENT ORDER EIGENVALUES AND EIGENVECTORS; 
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FOR I;=l UNTIL N DO 

BEGIN K: = I ; P: =D (I ) ; 

FOR J;=I+1 UNTIL N DO 
IF D( J ) < P THEN 

BEGIN K:=J;P:=D(J) ;END; 

IF K 1= I THEN 
BEGIN 

D( K) :=D( I ) ;D( I ) :=P; 

FOR J;=l UNTIL N DO 

BEGIN P:=Z(J,I) ;Z(J,I):=Z(J,K);Z(J,K):=P;END 

end:end; 

FIN; 

END TQL2; 

START:READ(N,NPX,NPY,NX,NY,HX,HY,EPS) ; 

BEGIN 

COMMENT PROGRAM STARTS HERE; 

REAL ARRAY D1 , D2 , El , E2 ( 1 : : N ) ; 

REAL ARRAY D , R ( 1 : : N, 1 : : N ) : 

REAL ARRAY U ( 1 : : N, r>: ; NPX , 0 : : NPY ) ; 

REAL ARRAY E t A , B , P ( 1 : : N, 1 : : N ) ; 

REAL ARRAY X(U::NPX) ; 

REAL ARRAY Y(0: : NPY) ; 

REAL ARRAY C , C 1 , C2 ( 1 : : N , 0 : i N PX ,0 : : NPY ) ; 

REAL ARRAY Z1 , Z2 ( 1 : : N , 1 : : N ) ; 

PROCEDURE FCN; 

COMMENT COMPUTES FUNCTIONS AND MATRICES A AND B; 

BEGIN 

X(0) :=Y(0) ;=0.0; 

FOR l;=l UNTIL NPX DO X(I):=I*HX; 

FOR l; = l UNTIL NPY DO Y ( I ) : = I=!=HY ; 

FOR I ;=1 UNTIL N DO 

FOR J;=l UNTIL N DO A ( I , J ) := B ( I , J ) : =0. 0 ; 

FOR I ;=0 UNTIL NPX DO 
FOR J:=0 UNTIL NPY DO 
BEGIN 

C(1,I,J): = (X(I)-X(I )**2 )^(Y( J)-Y( J )=!'*2 ) ; 

C(2 , 1 , J) ; = (0. 5-X( I ) )*C (1 , I , J) ; 

C(3,I , J) : = ( 0. 5-Y( J ) )*C( 1, I , J) ; 

Cl(ltI»J):=(l.-2.*X(I ) )*( Y( J)-Y( J )**2 ) ; 

C1(2»I,J); = (Y(J)-Y(J )=i'’t^2 0.5-3. ^X( I ) +3 .*X ( I ) ** 2 ) ; 

Cl(3,I,J); = (l.-2.-X(I) )*0.5-Y( J ) )*(Y( J )-Y( J )=f'=i=2) ; 

C2( 1, I , J) ; = ( l.-2.=NY( J ) )*( X( I )-X( I )*=i'2 ) : 

C2( 2,1 , J) ; = (0.5-X( I ) )^ (X( I )-X( I )**2)*( l.-2.*Y( J ) ) ; 
C2(3,I , J) ; = (X( I )-X( I )«#2)*(u. 5-3.*Y(J )+3.«Y( J )*=!‘2) : 
END; 

FOR K:=l UNTIL N DO 
FOR l;=i until n do 

BEGIN 

FOR I :=1 UNTIL NPX DO 
FOR J:=l UNTIL NPY DO 
BEGIN 

A(K,L) :=A(K,L )+Cl (K, I ,J)*d ( L , I , J ) X^HX^^HY 
+C2(K,I , J) *C2 ( L, I , J ) -HX*HY; 

8(K,L) ;=B(K,L)+C(K,I, J)*C(L, I,J)#HX*HY: 

END; 

END; 

END FCN; 

PROCEDURE WRITER; 

COMMENT WRITES ALL DATA; 

BEGIN 

FOR I;=l UNTIL N DO 

FOR J;=l UNTIL I DO 

BEGIN D0G;=Z1( I , J) ;Z1(I,J);=Z1(J,I);Z1(J,I) ;=COG; END; 
FOR L;=l UNTIL N DO 

FOR I; = l UNTIL N DO D ( L , I ) ; = Z2 ( I , L ) ; 

FOR L;=l UNTIL N DO 

FOR l;=l UNTIL N DO 

BEGIN D0G;=0.0: 

FOR J;=l UNTIL N DO DOG : = DOG+ Z1 ( I , J ) * D ( L , J ) ; 

R( I ,L) ;=DOG; 

END; 

FOR I ;=1 UNTI L N DO 
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BEGIN 

FOR J:=0 UNTIL NPX DO 
FOR L:=0 until NPY DO 

BEGIN U( I ,JtL) :=0.0; 

FOR K:=l UNTIL N DO U ( I , J » L ) : =U ( I , J , L ) +R ( I , K )* 

C(K, J,L) : 

end;end; 

FOR K:=l UNTIL N DO 
BEGIN IOCONTROL( 3) ; 

WRITE!" "); WRITE!" "); WRITE!" "); 

WRITE!" COMPUTER OUTPUT 5 . ! RAYL E I GH ) " ) ; 

WRITE!" ");WRITE!" ");WRITE!" "); 

WRITE!" EIGENFUNCTION=",K," OMEGA= " , SORT ! D2 ! K ) ) ) 

WRITE!" "); WRITE!" ") -.WRITE!" 

WRITE!" THE VALUES OF C!I)"); 

WRI TE!" ") ; WRITE! " ") ; 

FOR L:=1 until N DO 

BEGIN INTFIELDSIZE: = 1 ; 

WRITE!" C",L ,"=",R!K,L) ) ; 

WRITE!" ");WRITE!" " ) : END ; WR IT E ! " "); 

WRITE!" X Y U!X,Y)"); 

WRI TE!" ") :WRITE! " ") ; 

FOR J:=0 STEP NX UNTIL NPX DO 
FOR L:=0 STEP NY UNTIL NPY DO 
BEGIN INTFIELDSIZE:=2; 

WRITE!" ",J DIV NPX,".",J REM NPX»" ", 

L DIV NPY,".",L REM NP Y , U ! K , J , L ) ) ; 

WRITE!" ") ; 

END ; end: 

GO TO Q; 

END WRITER; 

FCN: 

TRED2!N,D1,E1 ,B,Z1,EPS) ; 

TQL2!N,DL,E1,Z1,T) : 

IF T THEN BEGIN WRIT E ! "F A I L= 1" ) ; GO TO Q; END ELSE 
BEGIN 

FOR I:=l UNTIL N DO D1 ! I ) : =SQRT ! D1 ! I ) ) : 

FOR I:=l UNTIL N DO 

FOR J:=l UNTIL N DO Z 1 ! J , I ) : =Z1 ! J , I ) / D1 ! I ) ; 

END ; 

FOR l:=l UNTIL N DO 
FOR K;=l UNTIL N DO 
BEGIN P!I,K):=n.O: 

FOR J:=l UNTIL N DO P ! I , K ) : =P ! I , K ) +A ! I , J ) *Z 1 ! J , K ) ; 

END; 

FOR I : =1 UNTI L N DO 
FOR J:=l UNTIL I DO 

BEGIN DOG:=Zl !I,J);Z1!I,J):=Z1!J,I);Z1!J,I) :=DCG; END; 

FOR I:=l UNTIL N DO 
FOR K;=l UNTIL N DO 
BEGIN E! I ,K) :=0.0 ; 

FOR J: = l UNTIL N DO E ! I , K ) : =E ! I , K) +Z 1 ! I , J ) =«= P ! J , K ) ; 

END ; 

TRED2!N,D2,E2,E,Z2,EPS) ; 

TQL2!N,D2,E2,Z2,T) : 

IF T THEN BEGIN WR I T E ! "F A I L=2 " ) ; GO TO Q;END; 

WRITER; 

END ; 

Q:G0 TO START; 

END, 
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COMPUTER PROGRAM 6. 



DYNAMIC PROGRAMMING APPLIED TO 

(1) V^U = 2(x^+y=^-x-y) EQ. (1.5) 

(2) V“U = X^U EQ. (1.3) 

BEGIN 

CCMMENT SOLVES PARTIAL DIFFERENTIAL EQUATIONS AND 
EIGENVALUE PROBLEMS ON A RECTANGLE BY DYNAMIC 
PROGRAMMING AND STODOLA'S METHOD. 



FORMAL PARAMETERS: 




QM 


OMEGA 




0M2 


OMEGA SQUARED 




N 


NUMBER OF X POSITIONS 




M 


NUMBER OF Y POSITIONS 




K1 


COEFF IF NEEDED 




VI 


ORDER OF THE OPERATOR 






V1=0 FIRST ORDER EIGENVALUE 


PRB 




Vl=2 SECOND ORDER EIGENVALUE 


PRB 


T1 


TYPE OF PROGRAM RUN 
T1=0 FIRST ORDER 
Tl=2 SECOND ORDER 
Tl=50 EIGENVALUE 




P( I ) 


SECOND NORMAL DERIVATIVE ON BDRY 


El 


FUNCTION TO SATISFY 




U 


UNKNOWN FUNCTION 




ei 


NUMBER OF EIGENVALUES; 





REAL K£,PE; 

REAL E ,ERA,OM,OM2 ,0M02 ,H ,K1 ; 

REAL CATl, CAT2.Cl,C2iC3.DOGl,DQG2; 
REAL ZNORMtZS; 

INTEGER N,M,N1,S,T1, VI ,Ml; 

INTEGER 01. Q2; 

LOGICAL TOGGLE; 

REAL ZNORMl; 

X:RSAD(N,M,H,K1,T1,Q1,V1); 

TOGGLE :=TRUE; 

0 2 : = 1 ; 

CAT1:=CAT2:=5000.0; 

Ml:=2*M-2; 

0M2:=1000.O; 

N1 : =0 ; 

ERA :=0.03; 

ERA:=0.002; 

E:=0.000001; 

BEGIN 

REAL ARRAY Xl( O: :N) ; 

REAL ARRAY Y1 (0: :M) ; 

REAL ARRAY U5 ( 0: :N ,0 : : M) ; 

REAL ARRAY U 1 . U2 ( 0 : : N . 0: : M ) ; 

REAL ARRAY U3 ( C : :N ,0 : : M ) ; 

REAL ARRAY C I ,B1 ,P1( 0: :N ,0: : M) ; 

REAL ARRAY Z(0 : : M+4, 0: :M+4) ; 

REAL ARRAY El ( 0 : : N ,0 : : M ) ; 

REAL ARRAY 1 1 , Q( 1 : : M-1 , 1 : : M- 1) ; 

REAL ARRAY U , 3 , R ( 0 : : N , 0: : M ) ; 

REAL ARRAY A ( 1 : : N , 1 : : M , 1 : : M ) ; 

REAL ARRAY D( 1 : :N. 0: :M , 0: :M1 ) ; 

REAL ARRAY CAT(0: :M) ; 

REAL ARRAY C ( 0 : : N ,0: : M ) ; 

REAL ARRAY PI, P3 (0 : : M ) ; 

REAL ARRAY P2,P4(0::N); 

FOR I:=0 UNTIL M DO READOMt U(0, 1 ) ) *, 

FOR T;=0 UNTIL N DO R5 ADON ( U ( 1 , 0 I ) ; 

FOR I:=0 UNTTL M DO READON( U( N, I ) ) ; 
FOR I:=0 UNTIL N DO R£ ADOM( U( I , M ) ) ; 

FOR 1:=0 UNTIL N DO XI ( I ) := I ; 

FOR J:=0 UNTIL M 00 Y1 ( J ) ; = J-*=H ; 

IF (T1>0) AMO (Vl=2) THEN 
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=0 until m do 
= 0 . 0 ; 

=0 UNTIL M DO 
=-Pl( I ) ; 

=0 UNTIL N DO 

“ 0 • 0 5 

=0*UNTIL N DO 
=-P2 ( I ) ; 



e=GiN 
FOR I 
Did) 

FOR I 
P3( I ) 

FOP I 
P2d ) 

FOR I 
P4 { I ) 

CNjn ; 

F0P“r:=0 UNTIL N DO 

FOR J:=0 UNTIL M DQ U1 ( I , J ) : =U2 ( I, J ) : =U3( I , J ) : =U5{ I , J ) : =0. 0 

S;=N; 

BEGIN 

PROCEDURE WRITER; 

BEGIN WR ITE( "SINGULAR") ; 

GO TO W; 

END WRITER; 

PPCCEDURS WPITl; 

BEGIN 

I0CCNTPGL(3) ; 

WRITEC ");WRITE(" ");WRITS(" ") ; 

WRITEC* COMPUTER OUTPUT 9. < DYNAMIC )") ; 

WRIT“(" ");WRITE(" ");WRITE(" "); 

IMTFIELDSIZE:=1; 

WRITEC EIGENFUNCTI0N = ",02-1 ) ;WRITE(" "); 

INTFIELDSI ZE : =2 ; 

WRITEC PATH = ",N1." 

WRITE (" ") ;WRITE(" " ) ; 



OMHGA=".OM) ; 

Y ", 



WRI TE ( " 

" U ( X , Y ) " ) ; 

WRITE(" "); WRITE (" " ) ; 

FOR l:=l STEP 2 UNTIL N-1 DO 
FOR J:=l STEP 2 UNTIL M- 1 DO 
BEGIN INTFISLDSIZS : =2 ; 

WRI TEC ",I DIV N,".".I REM N." 

J DIV M,".", J REM M,U( I,J ) ) ; 

WRI TE ( " ") ; 

END; 

IF Q2 = 01+l THEN GO TO W ELSE STARP ; 

END WPITl; 

PROCEDURE STARP; 

BEGIN 

IF Tl=50 THEN 
BEGIN 

IF 02=1 THEN 

BEGIN Q2:=2; 

PQR l;=l UNTIL N-1 DO 

DOR J:=l UNTIL M- 1 00 

U( I , J ) :=X1( I )=i^Yl( J)t. ( 1. -XKI ) )* (l.-Yl ( J ) ) ; 
GO TO Z5; 

END; 

IF 02=2 THEN 
BEGIN 

FOR I;=l UNTIL N-1 DO 

FOR J:=l UNTIL M- 1 DO 

BEGIN 

U2( I , J) ;=U( I, J) ; 

U5(I , J ) : = 0.0; 



U(I , J ) : =X1 ( I )'^(1.-X1( I ) ):^‘(. 5-XK I ) )*Y1( J )^( l.-YK J ) ) 
END; 

CATl: =KE; 02 : = 3;Nl:=l;0M2: =5000.0; 

GO TO Z5; 

END; 

IF 02=3 THEN 
BEGIN 

FOR I :=1 UNTI L N-1 DO 
FOR j:=l UNTIL M-1 DO 
BEGIN 

U3( I, J) :=U(I , J) ; 

U5(I,J):=0.0; 

U(I , J ) : =Y1 ( J) Y ( 1. -Y1 ( J ) )*( .5-Yl ( J ) )4=X 1 ( I ) ( 1 . -X 1 ( I ) ) 
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END ; 

CAT 2: =KE;02: = A; 0M2: = 5000.0; 

GO TO Z5; 

END; END; 

15: 

IF Q2>2 THEN GO TO Z6; 

COMMENT OUT FUNCTION TO SAT. HERE; 

IF Tl-=50 THEN 
BEGIN 

FOR I :=0 UNTIL N DO 
FOR J:=0 UNTIL M DO 
BEGIN 

El( I t J ) : = 2.*( XI (I )=^*2+Yl(J)*f2-Xl(I )-Yl( J) ) 
C(I ,J) :=2.0*H=^H*El(It J ) ; 

END; END; 

FOR l:=l UNTIL M-1 DO 
FOR J:=l UNTIL M-1 DO 
BEGIN 

IF I=J THEN I1(I.J): = 1.0 ELSE 1 1 ( I . J ) : =0. 0; 

Q(I ,J) :=0.0; 

IF I=J THEM Q(I ,J) :=K1»3.0; 

IF ABS(I-J)=1 THEN Q(I,J):=-K1; 

END; 

FOR j:=l UNTIL N DO 

FOR I;=l UNTIL M-1 DO 

BEGIN 

R( J .1 ) : =0. 0; 

IF 1=1 THEN R(J.I );=-2.0*Kl*U( J,0) ; 

IF I = M-1 THEN R( J ,I ):=-2.0^-Kl=«^U( J,M) ; 

END; 

FOR I:=l UNTIL M-1 DO 

FOR J:=l UNTIL M-1 DO 

A(M, I , J) : = I 1( I , J) ; 

FOR I: =1 until M-1 DO 
B(N,I ) :=-2.0*U(N,I ); 

SWITCH; 

Z6:!F Q2=2 THEN SWITCH ELSE OOP; 

END STARR; 

PROCEDURE OOP; 

COMMENT NORMALIZES U(X,Y); 

BEGIN 

ZN0RM:=0.0 ; 

FOR l:=0 UNTIL N DO 
FOR J:=0 UNTIL M DO 
BEGIN 

ZS: =ABS(U( I , J ) ) ; 

IF ZS>ZN0RM then ZN0RM:=ZS; 

END; 

FOR I:=0 UNTIL N DO 

FOP J:=0 UNTIL M DO U( I , J) : = U( I , J) / ZNORM; 
0MC2:=0M2; 

0M2:=1.0/ZM0PM; 

OM: =S0RT(0M2) ; 

IF Nl>30 THEM 
BEGIN 

WRITE (’'TOO MANY STEPS”); 

GO TO W; 

END; 

ZNORMl:=0.0; 

FOR I;=l UNTIL N-1 DO 
FOR J:=l UNTIL M-1 DO 
BEGIN 

ZS: =ABS(U( I , J )-U5 ( 1 , J ) ) ; 

IF ZS>ZN0RM1 THEN ZN0RM1:=ZS; 

END; 

IF znormkera then 
BEGIN KE;=0.0; 
f=OR I: = l UNTIL M DO 

FOR J;=l UNTIL DO KE : =KE+ (U ( I , J )*H ) «^2 ; 

WRITl; 

END; 

Nl;=Nl+l; 
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FOR I: =0 UNTIL N DO 
FOR J:=0 UNTIL M DO 
BEGIN 

U5( I , J) :=U( I , J) ; 

IF V1=0 THEN S1(I,J ): = -U( I,J ) ; 

IF Vl=2 THEN 51 (I , J) : =U( I , J) ; 

C(I»J):=2.*H*HAE1( I »J) ; 

END; 

FIGURE ; 

END DOP; 

PROCEDURE FIGURE; 

BEGIN 

comment solves for U AMD B; 

IF (T1=0) OR ((Tl=50) AND (V1=C)) THEN 
BEGIN 

FOR L:=M STEP -1 UNTIL 2 DO 
FOR I:=l UNTIL M-1 00 
BEGIN CAT( I ) :=0.0; 

FOR J:=l UNTIL M-1 DO 

CAT (I ) :=CAT(I )+D( L,I, J) *( B(L, J )+R(L-l, J )+C(L-l, J ) ) ; 
B(L-1,I ) :=CAT(I ) ; 

END; 

FOR I : = 1 UNTIL N-1 DO 
BEGIN 

FOR J:=l UNTIL M-1 DO 
BEGIN CAT( J) :=0.0; 

FOR L :=1 UNTIL M-1 DO 

CAT(J):=CAT(J)+D(I+1, J,L) 

’<'{U(I-1 ,U-(B(I + 1,L) + R(I,L)+C(I,L) )/2.0); 

U(I»J ):=CAT(J ) ; 

END; END; END ELSE 
BEGIN 

IF toggle then 
begin 

U1(0,0) :=-Pl( 0)-P2(0) ; 

U1(N,0):=P3(0)-P2(N); 

UKOtM) :=-Pl( M)+P4(0) ; 

UK N»M):=P4(M)+P3(M) ; 

FOR I:=l UNTIL M-1 DO 
BEGIN 

U1(0» I ) :=-Pl( I )+{ U(0. I+1)-2.^U(0.I )+U(0, I-l) ) /(H*H) ; 
U1 (N, I ) : = P3 ( I ) + (U( N, I+l)-2.-^'U(N, I )+U(N, I-l ) )/(H-’^H) ; 
END; 

POR I: = l UNTIL N-1 DO 
BEGIN 

UK I ,0) :=-P2(I ) + ( U(I + 1,0) -2.^'Ud ,vD)+U( I-l »0 ) )/{H^'H) ; 

U1(I,M); = P4(I )+(U(I+l,M )-2.=5^U( I,M)+U( I-1,M) )/(H^H) ; 

END; 

FOP I: = l UNTIL M-1 DO B1 ( N, I ) :=-2.*Ul ( N .1) ; 

FOR J:=l UNTIL N DO 
FOR I: = l UNTIL M-1 DO 
BEGIN RK J.I ) : = n.o; 

IF 1 = 1 THEN RKJ, I ): = -2.*Ul(J,0) ; 

IF I = M-l then R1( J,I) :=-2.^^U1 { J,M); 

END; 

TOGGLE :=FALSS; 

END; 

POR L:=N STEP -1 UNTIL 2 DO 
FOR I:=l UNTIL M-1 DO 
BEGIN CAT{ I ) : =0.0; 

FOR J:=l UNTIL M-1 DO 

CATC ) :=CAT(I > + 0(L, K (L, J)+R1(L-1, J )+C(L-l, J ) ) 

BKL-1,1) :=CAT( I ) ; 

END; 

FOR I:=l UNTIL M-1 DO 
BEGIN 

FOR J:=l UNTIL M-1 DO 
BEGIN CAT( J) :=0.0; 

FOR L:=l until M-1 DO 
C AT { J ) : =C AT { J ) -i-D ( I + 1 , J , L ) 

*( UK I-1,L) - ( RKI + l ,L)+R1 (I .L)+C( I .L) )/2. ) ; 

UK I, J ):=CAT(J ) ; 
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?MD;^ND; 

FOP I;=0 UNTIL N DO 

FOR J:=0 UNTIL M DO C I ( I , J ) : = 2 . *H*H*U1 ( I , J ) ; 

FOR L:=N ST^P -1 UNTIL 2 DO 
FOP. I: = l UNTIL M- 1 DO 
BHGIN CAT(I):=0.0; 

= OP J: = l UNTIL M-1 DO 

CAT (I ) :=CAT( I)+D(L»I.J)^(3(L,J)+R(L-1,J)+CI(L-1,J)); 

B(L-1 ,I ):=CAT(I ); 
cND; 

FOR I:=l UNTIL N- 1 DO 
BEGIN 

FOR J; = l UNTIL M-1 DO 
BEGIN CAT (J ) :=0.0; 

FOR L:=l UNTIL M-1 DO 
CAT(J ) ;=CAT( J )+D( I+l, J»L) 
*(U(I-l,L)-(B(I+l,L)+R(I,L)+CI(I,L))/2.); 

U( I t J) :=CAT( J) ; 

ENC;fnd;fnd; 

IF C2>2 then 
BEGIN 

Cl : = COGl : = DOG2*. =0.0; 

'=0R I: = l UNTIL N-1 DO 

FOR J: = l UNTIL M-1 DO C 1 : =C1 + U ( I , J ) 

FOR l:=l UNTIL N-1 00 

FOR J:=l UNTIL M-1 DO U( I , J ) : =U ( I , J ) -Cl ; 

FOR I:=l UNTIL N-1 DO 
FOR J:=l UNTIL M-1 DO 
BEGIN 

DOGl:=DOGl+U( I . J) vU2( I . J)=^H¥H; 

D0G2 : =D0G2 +U ( I , J ) =U3 ( I , J ) ; 

END; 

C2:=D0G1/CAT1 ;C3:=D0G2/CAT2; 

FOP, I;=l UNTIL N-1 DO 
FOR J: = l UNTIL M-1 DO 

U( I ,J ) : =U( I , J )-C2^U2( I ,J )-C3*U3( I , J ) ; 

END ; 

IF Tl=50 THEN DOP; 

I0CCNTR0L(3) ; 

WPITE(” 'M;WRITE(" •M;WRITE(" 

WPITEC’ computer output 6. (DYNAMIC)") 

WRITEC ");WRITE(" ");WRITE(" "); 

WR I TH ( " X Y " , 

•• U ( X , Y ) " ) ; 

WRI TE ( " ") ; WOI T£ (" " ) ; 

FOP I:=l STEP 2 UNTIL N-1 DO 
'^OR J:=l STEP 2 UNTIL M-1 CO 
BEGIN INTFIELDSIZ£:=2 ; 

WRITE!" "»I DIV N.",".I REM N." "♦ 

J DIV M, , J REM M,LI( I , J ) ) ; 

WRITE!" "); 

END; 

GO TO w; 

END FIGURE; 
procedure SWITCH; 

BEGIN 

INTEGER L,K; 

REAL DOM,AMUL,T; 

FOR L:=N STEP -1 UNTIL 2 DO 
BEGIN 
S; = S — 1 ' 

FOR I: = l UNTIL M-1 DO 
FOR J:=l UNTIL M-1 DO 

D!L, I , J ) :=Q! I, J )+A!L, !, J ) ; 

FOR I : = 1 UNTIL M-1 DO 
FOR J:=l UNTIL M-1 DO 
D!L,I,J + M-1 ):=I1!I,J); 

= 0R K; = l UNTIL M-1 DO 
PFGIN 

IF K=M-1 then go to P; 

FOR I:=K+1 UNTIL M-1 DO 
BEGIN 
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IF ABS(D(L,K,K ) )<ABS( D(L, I,K ) ) THEN 
BEGIN 

FOR j;=l UNTIL 2*M-2 DO 
BEGIN 

T;=D( L,I , J) ;D(L,I ,J) :=D(L,K, J); 
D(L.K,J):=T; 

END;END;END; 

P:IF ABS(D(L,K,K) )<E THEN WRITER ELSE 
BEGIN 

DOM:=D(L,K,K) ; 

FOR J:=l UNTIL 2*M-2 DO D ( L . K , J ) : =D ( L . K , 
FOR I:=l UNTIL M-1 DO 
BEGIN 

AMUL:=D(L,I ,K) ; 

IF I=K THEN ELSE 
BEGIN 

FOR J:=l UNTIL 2*M-2 DO 
BEGIN 

D(L,I,J):=D(L,I,J)-AMUL*D(L, 

end;END;snd;END;ENO; 

FOR I:=l UNTIL M-1 DO 
FOR J: = l UNTIL M-1 DO 
D(L. I.J ) : = D(L» I » J+M-1) ; 

FOR !:=1 UNTIL M-1 DO 
FOR J:=l UNTIL M-1 DO 

A(L-1. I.J ):=I1( I,J)-0(L,I,J) ; 

END; 

IF Tl=50 THEN DOP ELSE FIGURE; 

END SWITCH; 

STARP; 

END; 

END; 

W:G0 to X; 

Z:END. 



J )/DOM 



K,J ) ; 
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COMPUTER PROGRAM 7- 



DYNAMIC PROGRAMMING APPLIED TO 



(1) V^U = -X^U EQ. (1.1) 



(2) V**U = 8.0 EQ. (1.7) 

BEGIN 

COMWHNT SOLVES PARTIAL D I FPERENT I AL EQUATIONS AND 
EIGENVALUE PROBLEMS ON A RECTANGLE BY DYNAMIC 
PROGRAMMING AND STODOLA’S METHOD. 

FORMAL PARAMETERS: 

CM OMEGA 

CM2 OMEGA SQUARED 

N NUMBER OF X POSITIONS 

M number of Y POSITIONS 

K1 COEPF IF NEEDED 

VI ORDER OF the OPERATOR 

V1=0 FIRST ORDER EIGENVALUE PPB 
Vl=2 SECOND ORDER EIGENVALUE PRB 
T1 TYPE OF PROGRAM RUN 

T1=0 FIRST ORDER 
Tl=2 SECOND ORDER 
Tl=50 EIGENVALUE 

P(I) SECOND NORMAL DERIVATIVE ON 6DRY 
El FUNCTION TO SATISFY 

U UNKNOWN FUNCTION 

Q1 number of EIGENVALUES; 

REAL KEiPE; 

REAL S,SRA,0M,0M2,0M02,H,K1 ; 

REAL CATl ,CAT2 ,C1 ,C2 ,C3, DOGl ,DOG2; 

REAL ZNORM.ZS; 

INTEGER N,MtN1tS,T1,V1,M1; 

INTEGER Ql,02; 

LOGICAL TOGGLE; 

REAL ZNORMl; 

X:PEA0(N,M,H,K1,T1 ,Q1 ,V1 ) ; 

TOGGL5:=TRU5; 

02 : = 1 ; 

CATl:=CAT2:=5000.0; 

Ml: =2’!- ^^-2; 

0M2 :=1 000.0; 

Nl:=0; 

ERA:=0.03; 

ERA :=C .002; 

E:=0. 000001 ; 

BEGIN 

REAL ARRAY XI (0: :N); 

REAL ARRAY Y1 ( 0: :M) ; 

REAL ARRAY U5 ( 0: : N .0 : : M ) *, 

REAL ARRAY U1 » U2 (0 : : N ,0 : : M ) ; 

REAL ARRAY U3 ( 0: : N ,0 : : : 

REAL ARRAY Cl , B1 , R1 ( 0: :N , 0: : M ) ; 

REAL ARRAY Z ( 0 : : N+4 , 0 : : M+A ) ; 

REAL ARRAY El ( 0: . 0 : : M) ; 

REAL ARRAY I 1 , Q{ 1 : :M-1 , 1 : :M- 1 ) ; 

REAL ARRAY U , 3 ,R ( 0 : : N ,0 : : w) ; 

REAL ARRAY A (1 : : N » 1: : M , 1 : : M ) ; 

R^AL ARRAY D ( 1 : : M tC : : M ,0 : : Ml ) ; 

REAL ARRAY CAT(0: :M) ; 

REAL ARRAY C ( 0 : : N, 0: : M ) ; 

REAL AR^AY pi ,P3(0:: M) ; 

REAL ARRAY P2.P4(0: :N) ; 



FOR I 
FOP I 
FOR ! 
FOP I 
FOR I 
FOR J 



=0 UNTIL M DO READ0M(U(0, I ) ) ; 

=0 UNTIL N DO READON(U(I tO) ) ; 

=0 UNTIL M DO '^HAD0N( U( N, I ) ) ; 

=0 UNTIL N DO R- ADONUJ ( I ,M ) ) ; 

=0 UNTI L N DO XI (I ) :=H’H; 

=0 UNTIL M DO Y1 ( J ) :=J^H; 



IF (T1>0) AND (Vl=2) THEN 
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B^GIN 
FOR I 
PI ( I ) 
FOP I 
P3 n ) 
FOP I 
P2( I ) 
FOP 1 
P4(!) 



=0 UNTIL M DO 

=2,^ { Y1 (I )-Yi ( I )=*=*2 ) ; 

=0 UNTIL M DO 

=-Pl n ) ; 

=0 UNTIL N DO 
= 2.^ ( Xl( I )-Xl( I )+*2) ; 
=0 UNTIL N DO 
=-P2( I ) ; 



END; 

FOR I:=0 UNTIL N DO 

FOR J:=0 UNTIL M DO U1 ( I t J) : =U2 ( I , J) : =U3 ( I , J ) : =U5 ( I . J ) : =0 .0 
S: = N; 

BEGIN 

PP.CCEDUP5 WP.ITFR; 

BEGIN WRITSC'SINGULAR" ) ; 

GO TO W; 

END WR ITHR ; 

PROCEDURE WRITl; 

BEGIN 

IOCONTROL( 3) ; 

WRITHC »);V!RITE(" '');WRITEC' "); 

WRITSC COMPUTER OUTPUT 7 . ( DYNAMI C ) " ) ; 

WRITF(» ••); WRITS!" ");WRITS(" •• ) ; 

INTFIHLCSIZE:=l; 

WRITE!" £IGENFUNCTI0N=",02-1 ) ;WRITH!" 

INTPIELDSIZE:=2; 

WRIT”!" PATH=",N1," OMEGA="»OM); 

WRITE! " "); WRITE ! " ") ; 

WRITS!" X Y "t 



" U!X,Y)"); 

WRITE! " "); WRITS! " ") ; 

FOR I;=l STEP 2 UNTIL N-1 CO 
FOR J:=l STEP 2 UNTIL ^<-1 DO 
begin INTFrELDSIZE;=2 ; 

WRITE!" ",I DIV N,".",I RSM N," "t 

J DIV RSM M,U!I ,J) ); 

WRITS!" " ) ; 

END; 

IF 02=01+1 THEN GO TO W ELSE STARP; 

END WRITl; 

PROCEDURE STARP; 

BEGIN 

IF Tl=50 THEN 
BSGIN 

IF 02=1 THEN 

BSGIN 02: =2; 

FOR I:=l UNTIL N-1 DO 

FOR J:=l UNTIL M- 1 DO 

U!I ,J) :=X1!I)*Y1! J)-^!1.-X1 !I) )*!1,-YUJ) ) ; 

GO TO Z5; 

END; 

• IF 02=2 THEN 
BSGIN 

FOR l:=l UNTIL N-1 DO 

FOR J;=l UNTIL V-l DO 

BSGIN 

U2!I,J):=U! I,J); 

U5! I , J) ;=0.0; 

U!I,J): = Xl!I)F:(l.-Xim)^!.5-Xl!I ) ) *Y 1 ! J ) i' ! l.-Yl !J ) ) 
END; 

CAT1:=KE;Q2:=3;N1 :=l; 0M2:=5000.0; 

GO TO Z5; 

END; 

IF 02=3 THEN 
BEGIN 

FOR I :=1 UNTI L N-1 DO 
FOR J;=l UNTIL M-1 DO 
BEGIN 

U3! I , J) :=U! I , J ) ; 

U5! I, J ) ; = 0.0; 

U!I , J ) : =Y1 ! J)»M 1. -Y1 ! J ) )*! .5-Yl! J ) )4=X1! I )^! l.-Xl! I ) ) 
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END; 

CAT2:=KE ;02:=4;0M2:=5000. 0; 

GO TO Z5; 

END; END; 

Z5: 

IF Q2>2 THEN GO TO Z6; 

COMMENT PUT FUNCTION TO SAT. HERE; 

IF Tl-.= 50 THEN 
BEGIN 

FOR I :=0 UNTI L N DO 
FOR J:=0 UNTIL M DO 
BEGIN 

Eld f J) : = 8.0; 

C(I.J):=2.0*H*H*E1(I.J) ; 

END; END; 

FOR I :=1 UNTIL M-1 DO 
FOR J:=l UNTIL M-1 DO 
BEGIN 

IF I=J THEN I1(I,J):=1.0 ELSE I1(!,J):=0.0 
Q(I .J ):=0.0; 

IF I=J THEN Q( I ,J ) :=K1T3.0 ; 

IF ABS(!-J)=1 THEN g(I,J):=-Kl; 

END; 

FOR J:=l until N do 

FOR I:=l UNTIL M-1 DO 

BEGIN 

R ( J , I ) : =0 . 0 ; 

IF 1 = 1 THEN R( J ,I ) : =-2.0*Kl=;*U( J,0) ; 

IF I=M-1 then R( j , ! ): =-2.0»^KliU(J ,M) ; 

END; 

FOR I:=l UNTIL M-1 DO 

FOR J:=l UNTIL M-1 DO 

A(NtI »J) :=i 1(1 »J); 

FOR I:=l UNTIL M-1 DO 
B(N,I ) ;=-2.0'-U(N,! ) ; 

SWI TCH ; 

Z6:IF Q2=2 THEM SWITCH ELSE OOP; 

END STARP; 

PROCEDURE OOP; 

COMMENT ^10RMALIZES U(X,Y); 

BEGIN 

ZNDRM: =0.0; 

FOR I:=0 UNTIL N DO 
FOR J:=0 UNTIL M DO 
BEGIN 

ZS:=ABS(U( I,J ) ); 

IF ZS>ZNORM THEN ZNORM;=ZS; 

END; 

FOR l:=0 UNTIL N DO 

FOR j:=0 UNTIL M DO U ( I , J ) : =U ( I , J ) / ZMORM ; 
OMC2:=OM2; 

GM2 :=1 .0/ZMDRM; 

OM:=SQf!T( 0^2) ; 

IF Nl>30 THHN 
BEGIN 

WRITE ("TOO MANY STEPS"); 

GO TO W; 

2ND; 

ZNORMl : = 0. 0; 

FOR I:=l UNTIL N-1 DO 
FOR J:=l UNTIL M-1 DO 
BEGIN 

ZS:=A3S(U( dJ )-U5( I,J )) ; 

IF ZS>ZN0RM1 THEN ZNORMl :=ZS; 

END; 

IF ZNORMKERA THEN 
BEGIN Kr:=0.0; 

FDR l:=l UNTIL N DO 

FOR J:=l UNTIL M DO K5 : = K £+ ( U ( I » J ) *H ) ; 

WRI Tl; 

END; 

N1 ; =N1+1 ; 
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FOR I:=0 UNTIL N DO 
FOR J:=0 UNTIL M DO 
BcGIN 

U5(I ,J) :=U(I , J) ; 

IF V1=0 THSN ?1( I .J ):=-U( I, J) ; 

IF Vl=2 TH'£N FI ( I , J ):=U( I , J ) ; 

C(I tJ) :=2.>^H*H«F1(I ,J) ; 

END; 

FIGURE; 

FND DOP; 

PROCEDURE FIGURE; 

BEGIN 

COMMENT SOLVES FOP U AND B; 

IF (T1=0) OR ((Tl=50) AND (V1=0)) THEN 
BEGIN 

FOR L:=N STEP -1 UNTIL 2 DO 
FOR l:=l UNTIL M- 1 DO 
BEGIN CATd ) :=O.C; 

FOP J: = l UNTIL M-1 DO 

CATd ) :=CAT ( I )+D(L, I» J )'^( B(L, J y+R (L-1 , J )+C( L-1, J) ) ; 
B(L-ld) :=CAT(I ); 

END; 

FOR l:=l UNTIL N-1 DO 
BEGIN 

FOR J;=l UNTIL M-1 DO 
BEGIN CAT(J):=0.0; 

FOR L: = l UNTIL M-1 DO 

CAT(J):=CAT( J)+D( I + l,JtL) 

*(Ud-l ,L)-('’(I+l,L)+R(IiL)+C(I,L) )/2.0); 

U( I , J ) :=CAT(J ) ; 

END;END;SND ELSE 
BFG I N 

IF toggle then 

BEGIN 

U1(0,0) :=-Pl(0)-P2(0) ; 

U1(N,0):=P3(0)-P2(M); 

U1 (0,M) :=-Pl ( M)+PA(G ) ; 

UK N,M) :=PA( N)+d3( M) ; 

FOR l:=l UNTIL M-1 DO 
BEGIN 

UK 0, 1 ) :=-Pl( I ) + ( U{0,I + 1) -2.* U(0,n+U(0 . I-l ) )/ ; 

U1 (N, I ) :=P3( I ) + (U (N, I + l)- 2 .^ni(N, I )+U( N , I-l ) ) ; 

END; 

FOP I: = l UNTIL N-1 DO 
BEGIN 

Uld ,C) :=-P2d ) + (Ud+l ,0 )-2 .<MJ( I ,0 )+U( 1-1,0 ) )/(H’^H) ; 
Uld ,M) : = P4( I ) + ( U( I + 1,M )-2.^ai( I ,M)+Ud -1 , M) ) /(HtH) ; 
END; 

FOR I; = l UNTIL M-1 DO B1 ( N , I ) ; =-2 .*U1 ( N , I ) ; 

FOR J:=l UNTIL N DO 
FCR I:=l UNTIL M-1 DO 
begin P1( J,I ) :=0. 0; 

IF 1 = 1 THEN Rl( J, I) :=-2.i-Ul( J ,0) ; 

IF I = M-1 them F.K J,I) :=-2.^U1 (J,M); 

END ; 

TOGGLE :=FALSS; 

END; 

FOR L:=N STEP -1 UNTIL 2 DO 
’^OR T:=l UNTIL M-1 DO 
BEGIN CAT(I):=C.O; 

FOR J:=l UNTIL M-1 DO 

CAT d ) :=CAT( I ) +0 ( L , I , J )* ( B1 ( L , J ) +R1 ( L- 1 , J ) +C ( L- 1 , J ) ) 
Bl( L-1 ,1) :=CAT( I ) ; 

END; 

FOR l:=l UNTIL N-1 DO 
BEGIN 

FOR J:=l UNTIL M-1 DO 
BEGIN CAT(J):=0.0; 

FOR L:=l LlfiTI L M-1 DO 
CAT ( J ) :=CAT ( J ) 40( I + l, J ,L ) 

t(Uld-l ,L)-( B1 (1+1 ,L)+R1 (I ,L)+C( I,L) )/2. ); 

U1 ( I, J ) :=CAT( J ) ; 
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srJD; rNO; 

FCP I : =0 UNTI L N 00 

POP J:=C UNTIL M DO C I ( I . J ) : = 2. ( I » J ) ; 

POR L:=M STEP -1 UNTIL 2 DO 
FOR I: = l UNTIL M- 1 DO 
BEGIN CAT(I):=0.0; 

FOR J;=l UNTIL K- 1 DO 

CaT(I):=CAT( I )+D(LtI ,J)’^(B(L.J) + P(L-1.J)+CI(L-1,J) ); 
P(L-1 ,I):=CaT(I ); 

END; 

FOR I; = l UNTIL N-1 DO 
B5GIN 

FOR J : = 1 UNTI L N-1 DO 
BEGIN CAT( J ): = 0.0; 

FOR L:=l UNTIL M-1 DO 
CAT( J) :=:aT( J)+D(I + ] , j,L) 

1=(U(I-1.L )-(8(I+l,L)+R( I,L)+CI(ItL) )/2. ) ; 

U( ! ,J ) :=CAT ( J ) ; 

2ND; END; END; 

IF C2>2 THEM 
BEGIN 

Cl:=DOGl:=DOG2:=C.O; 

FOR I:=l UNTIL N-1 DO 

FOR J: = l UNTIL M-1 DO Cl : =C1 +U ( I t J ; 

FOR l: = l UNTIL N-1 DO 

FOR J:=l UNTIL M-1 DO U ( I , J ) : =U ( I , U )-C 1 ; 

FOR I:=l UNTIL N-1 DO 
FOR J: = l UNTIL iM-1 DO 
BEGIN 

DOG1:=DOG1 + U(I , J) *^U2(I , J)*H4 H; 

DQG2:=D0G2+U( I , J )tU3( I,J )+H^H ; 

END; 

C2:=D0G1/CAT1 ;C3: =DGG 2/CAT2 ; 

FOR l:=l UNTIL N-1 DO 
FOP J;=l UNTIL M-1 DO 
U( I , J );=U( I. J )-C2TU2( I ,J)-C3*U3( I » J) ; 

END; 

IF Tl=50 THEN DOP; 
lOCCNTROL (3) ; 

WRITEC *');WPIT5C’ ");WRITE('« "); 

WRITEC COMPUTER OUTPUT 8, ( DYNAMIC )" ) 

VIRITEC »);WRITE(" » ) ;WRI T5 ( •* »); 

WRITE!” X Y ”, 

” U(X,Y)"); 

WPI Tc ( ” ”) ;WRITE (” ” ) ; 

FOR I:=l STEP 2 UNTIL N-1 00 
FCR J:=l STEP 2 UNTIL M-1 DO 
BEGIN INTFIEL0SIZE;=2 ; 

WRITE!” ”,I DIV N,”.”,I REM N»" 

J DIV M,”.”, J REM M,U!I,J) ) ; 

WRITE!” 

END; 

GO TC w; 

END FIGURE; 

PROCEDURE SWITCH; 

BEGIN 

INTEGER L,K; 

REAL DOM,AMUL.T; 

FCR L:=N STEP -1 UNTIL 2 DO 
BEGIN 
S ; = S- 1 * 

FOR I;=l UNTIL M-1 DO 

FOR J:=l UNTIL M-1 DO 

D !L, I, J ) : = Q ! I, J )+A! L , I , J ) ; 

FOR I;=l UNTIL M-1 DO 

FOR J:=l UNTIL M-1 DO 

D!L, I,J+M-1); = I1!T,J); 

FCR K:=1 until M-1 DO 
BEGIN 

IF K=M-l then GO TO P; 

FOR I:=K+1 UNTIL M-i OO 
BEGIN 
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IF ABS(0(L,K,K) )<ABS(D(L,I .K)) THEN 
BEGIN 

FOR J:=l UNTIL 2^M-2 DO 
BEGIN 

T:=C(L,ItJ);D(LtI»J):=D(LtKtJ); 

D ( L»Kt J) :=T; 

ENC;END;End; 

P:IF ABS(D(L,KtK) )<r THEN WRITER ELSE 
BEGIN 

DOM: = D(Lt K, K ) ; 

FOR J:=l UNTIL 2*M-2 DO D ( L , K t J ) : = D ( L , K , J )/ DOM 
FOR l:=l UNTIL M- 1 DO 
BEGIN 

AMUL: =D(L,I,K); 

IF I=K THEN ELSE 
BEGIN 

FOR J:=l UNTIL 2>^M-2 DO 
BEGIN 

D(L,ItJ):=D(L, It J)-AMUL=«‘0(LtK,J ) ; 
END;END;END;END;3ND; 

FOR I:=l UNTIL M-1 DO 
FOR J: = l UNTIL M-1 DO 
D(Lt I tJ ):=0(L» I, J+M-1) ; 

FOR I :=1 UNTI L M-1 DO 
FOP J: = l UNTIL M-1 DO 

A(L-1,I tJ )J =I1( It J)-D(L, ItJ ) ; 

END; 

IF Tl=50 THEN DOR ELSE FIGURE; 

END SWITCH; 

STARP; 

END ; 

END; 

W-.GQ TO X; 

Z: END. 
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fi 

a 



COMPUTER PROGRAM 8 



GALERKIN'S METHOD APPLIED TO 
V^U = PCx'^+y^'-x-y) EQ. (1.5) 



BEGIN 

COMMENT SOLVES PAP^T!AL DIFFERENTIAL EQUATIONS 0'= 

THE FORM L(U)=F BY GALERKIN'S METHOD ON ANY DOMAIN. 
FORMAL PARAMETERS; 

N NUMBER OF EQUATIONS 

M2 MAX. number OF X OQSITIONS 

M3 MAX. NUMBER OF Y POSITIONS 

N2 ARRAY OF LOWER Y POSITIONS 

N3 ARRAY OF UPPER Y POSITIONS 

F KNOWN FUNCTION 

U UNKNOWN FUNCTION 

D(I) ESTIMATING FUNCTIONS 
Dl(T) L(ESTIMATING FUNCTIONS) 

Ad) COSF. OP ESTIMATING FUNCTIONS; 

REAL H,E,DOG,CAT,DEV,AMUL,T,OOM; 

REAL Hi; 

REAL H2,H3; 

INTEGER M2,M3; 

INTEGER N,M,R; 

2 : =0. 00000 1 ; P : =0 ; DEV : =1 . 0 ; 

READ(N.M2.M3.H2.H3) ; 

HI; =H2*H3; 

BEG I N 

REAL ARRAY X(0.*;M2); 

REAL ARRAY Y<0;;m3); 

REAL ARRAY B( 1 ; ; M , 1 *• ; N+1 ) ; 

REAL ARRAY A»C(l; ;N) ; 

PEAL ARRAY F , U (0 ; M2 ,0 ; ; M3 ) ; 

REAL ARRAY D ,D1 ( 1 ; ;N ,0 ; ; m2 ,0 ; ; M3 ) ; 

INTEGER ARRAY N2 . N3( 0 ; ; M2 ) ; 



FOR I 
FOR I 
FOR I 
FOR I 
FOR I 
FOR 



0 UNTIL M2 DO X ( I ) ; = I-^H2 ; 

=0 UNTIL M3 DO Y(I);=I <H3; 

=0 UNTIL M2 DO R EADON ( N 2 ( ! ) ) ; 
=0 UNTIL M2 DO RFADON( N3 { I ) ) ; 
=0 UNTIL M2 DO 



=0 UNTIL M3 DO 
BEGIN 

P( I ,J) ;=0.0; 

U(I ,J);=0.0; 

FOR L;=l UNTIL N DO 

D(L,I,J) ;=D1(L.I ,J) ;=o.o; 

END; 

BEGIN CCMM'=NT SOLVE HERE; 

PROCEDURE WRITER; 

BEGIN WRITE( "S INGULAR" ) ; GO TO 0; 

END WRITER; 

PROCEDURi PERFORM; 

BEGIN COMMENT GET U(I,J) AND Ad); 

WRITEC ");WRITE(" '»);WRITE('' ") ; 

COMPUTER OUTPUT 10. <GALERKI M) ") 



WRITEC 
WRI TE( " ") ; WRITE ( " " ) ; 
WRI TE( " 

WRITEC ») ; 

FOR I ; =1 UNTI L N 00 
BEGIN 

Ad ) ; = B( I ,N+1 ) ; 
INTFIELDS! ZH;=l; 
WRITEC 
END; 

FOR J;=0 UNTIL M3-1 DO 
BEGIN 

FOR l;=0 UNTIL M2-1 DO 



VALUES OF C(I) ARE"); WRITEC "); 



C«,I ." = ".Ad ) ) ;WRITE c ") ; 
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B5GIN COG: =0.0; 

FOR L:=1 UMTIL N DO DOG: =OOG+A( L)*D ( L ,I . J) ; 

L( I ,J) :=DOG; 

5ND;0ND; 

WRITSC ");WRITE(" ”);WRITEC* ” ) *, 

WRIT5C' X Y '•» 

" U( X,Y)") ; 

WRITS (•' ") :WPITE(" " ) ; 

FOR I:=l ST'P 2 UNTIL M2-1 DO 
FOR J:=l STBP 2 UNTIL M3-1 DO 
BEGIN INTFI5L0SIZE: =2; 

WRITcC »SI DIV M2»»‘.'M REM M2, 

•• "tJ DIV M3,''.”,J REM M3,U(I,J)); 

WRI TE( " ") ; 

END; 

GO TO 0; 
cND PERFORM; 

PROCED'JPF switch; 

BEGIN COMMENT SOLVES FOR A(I) COEFF; 

INTEGER 

Ml:=N;Nl:=N+l ; 

FOP I:=l UNTIL Ml DO 3 ( I , Ml +1 ) : =C ( I ) ; 

FOR K:=l UMTIL Ml DO 

BEGIN IF K=M1 THEN GO TO P; 

FOR I:=K+l UNTIL Ml 00 

BEGIN IF APS( B( K, K) )<AES( B( I ,K)) THEN 
BEGIN R:=R+l; 

FOR J:=l UNTIL M DO 

BEGIN T:=B(T,J);B(I ,J):=B(K,J);B(K,J): = T;END 

EM D ; END ; 

P:IF ABS(5(K,K)XE THEN WRITER ELSE 

BEGIN DEV:=OS VyB( K,K) ;D0M:=B( K,K) ; 

FOR J:=l UMTIL Ml DO B ( K , J ) : =B ( K , J ) /DOM ; 

FOR I:=l UMTIL Ml DO 

BEGIN AMUL:=B( I ,K) ; 

IF I=K THEN ELSE 
BEGIN 

FOR J:=l UNTIL N1 00 

BEGIN B( I ,J ) :=B( I , J )-AMUL*B ( K , J ) ; 

END; END; END; END; END; 

PERFORM; 

£NC SWITCH; 

COMMENT GUESSED FUNCTIONS 0(1) HERE; 

FOP J:=0 UNTIL M2 DO 
FOR K:=N2(J) UNTIL N3(J) DO 
3EGI N 

F( J ,K ) :=2.-^( X( J)**2+Y(K)X:=!£2-X( J)-Y(K) ) ; 

0(1 , J,K) : = ( X( J )X*2-X( J )-*-*3)9( Y( K)-Y( K )t-*2) ; 

D( 2, J,K) : = ( X( J) -X( J )* ( Y ( K )iJ*2-Y ( K ) ^«3 ) ; 

D(3,J,K ): = ( X( J )*?2-X( J)^'^3)>^( Y( K) t5.2-Y( K)*^3) ; 

D1 (1 , J ,K) : = (2 .-6 ,^X ( J ) )«(Y(K ) -Y(K )^t2) + 

2.«( X( J ) + ^3-X( J )^--2) ; 

DK2, J,K) : = (2.-6.’^Y(K ) )>•( X(J)-X(J )*’V2) + 

2.Y(Y(K)X::*c3-Y(K)‘«-'.2) ; 

Dl(3,J,K):=(2.-6.yX(J) ) ? ( Y ( K) ’f-^2-Y ( K )^>i'3 ) + 

(2.-6.^Y(K ) )i«(X( J )**2-X( ; 

END; 

COMMENT SOLVE FOR INTEGRALS F^D AND LDX=D; 

FOR I : =1 UNTIL N DO 
BEGIN DOG:=C.O; 

FOR J:=l UNTIL M2 CO 

POR K:=N2(J)+1 UNTIL N3(J) DO 

DOG:=OOG+(F( J,K)>^D(I , J,K)*H1) ; 

C( I ): = DGG; 

FOR J:=l UMTIL N DO 
BEGIN CAT:=0. 0; 

FOR K:=l UNTIL M2 DO 

FOR L:=N2(K)+1 UNTIL N3(K) DO 

CAT:=CAT+(D1( J , K, L)*D ( I ,K ,L ) ^H1 ) ; 

B(I ,J ) :=CAT; 

END; 

IF I=N THEN SWITCH; 



99 



end; 

£ND;End; 

0:END. 
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